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PREFACE 


This  report  incorporates  the  work  done  in  a  number  of  different  efforts 
to  improve  the  Articulated  Total  Body  (AIB)  Model's  capability  to 
simulate  human  body  biomechanics  in  vanouE  dynamic  environments, 
especially  aircraft  ejection  with  windblast  exposure. 

The  majority  of  modifications  to  the  model  fall  into  six  categories: 

wi n d  f  orc  e  opt  ion 
joint  diift  correction 
edge  effect  option 
multi-axis  angular  displacement 
vehicle  motion  prescription 
slip  joint  option 
hy perel 1 1 psoi d  option 

These  improvements  have  been  combined  to  form  the  ATB-IV  version  on  the 
Armstrong  Aerospace  Medical  Research  Laboratory's  (AAMRL)  Concurrent 
computer  system  at  Wright  Patterson  Air  Force  Base.  AAMRL,  Systems 
Research  Laboratories  Inc..  J  &  J  Technologies  Inc.,  and  the  National 
Highway  Traffic  Safety  Administration  have  all  contributed  to  the 
technical  work  described  herein. 
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1.0  INTRODUCTION 


The  Articulated  Total  Body  (ATB)  Model  is  used  at  the  Armstrong 
Aerospace  Medical  Research  Laboratory  (AAMRL)  for  predicting  gross  human 
body  response  in  various  dynamic  environments,  especially  aircraft 
ejection  with  vindblast  exposure.  Aerodynamic  force  application  and  a 
harness  belt  capability  were  added  to  the  Crash  Victim  Simulation  (CVS) 
Program  (Ref.  1),  by  Calspan  Corporation  in  1975  for  AMRL  (Ref  2.),  and 
the  resulting  program  became  known  as  the  ATB  model.  In  1980,  Calspan 
made  a  number  of  modifications  to  the  ATB  model  combining  it  with  the 
then  current  3-D  Crash  Victim  Simulation  program  to  form  the  ATB-1I 
model  (Ref.  3).  Complete  documentation  of  the  program  through  the 
ATB-II  version  was  performed  by  Calspan  Corp.  (Ref.  4).  A  new  version, 
ATB-I1I,  wa6  generated  which  included  the  improvements  made  by  J  &  J 
Technologies  Inc  to  model  the  body  response  to  windblast  for  AMRL  (Ref. 
5). 


A  number  of  additional  efforts  have  been  made  to  improve  various  aspects 
of  the  ATB-III  model,  with  emphasis  on  its  capability  to  simulate 
aircraft  ejection  with  windblast  exposure  as  well  as  complex  automobile 
accidents . 

This  volume,  Modifications,  contains  a  description  of  the  major  changes 
made  to  create  ATB-1V  and  the  theory  U6ed  to  develop  them. 

Section  Two  of  this  volume  includes  a  new  wind  force  option  allowing 
segment  contact  ellipsoids  to  block  the  wind  as  well  as  other 
aerodynamic  force  improvements.  Corrections  to  prevent  angular  drift  in 
the  joints  are  described  in  Section  Three.  The  edge  effect  option  in 
Section  Four  ensures  that  a  contact  of  a  plane  with  an  ellipsoid  will 
not  be  ignored  and  that  a  smaller  force  will  be  applied  when  only  part 
of  the  contact  area  is  withiu  the  plane  boundaries.  Section  Five 
contains  an  improvement  allowing  the  prescription  of  multi-axis  angular 
displacements  to  describe  the  vehicle  motion.  A  new  option  allowing  a 
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joint  to  slice  along  an  axis  is  explained  in  Section  Six.  Section  Seven 
contains  a  new  hyperellipsoid  option.  A  summary  of  other  modifications 
that  form  the  ATB-IV  version  is  included  in  Section  Eight. 

These  changes  have  been  made  so  that  previous  input  decks  are  valid  with 
changes  required  only  in  the  H  cards.  The  updated  input  description 
outlining  any  changes  needed  and  describing  the  use  of  the  new  options 
is  described  in  Volume  2,  the  User's  Guide.  Sample  input  decks  using 
the  new  options  and  the  resulting  output  are  also  included  in  Volume  2, 
along  with  an  updated  list  of  numbered  stops.  Volume  3,  the 
Programmer's  Guide,  contains  the  listing  of  the  updated  ATB-IV  program. 
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2.0  AERODYNAMIC  FORCES 


The  aerodynamic  routines  added  to  the  model  in  1975  (Ref.  2)  have 
some  limitations  which  make  simulating  wind  forces  difficult.  For 
example  the  aerodynamic  pressure  is  prescribed  tabularly  as  a  function 
of  time,  but  this  requires  knowing  the  velocity  profile  of  the  seat 
before  a  simulation  is  made.  Since  the  3eat's  motion  may  depend  on  the 
wind  forces,  estimation  or  trial  and  error  has  to  be  used  in  defining 
the  aerodynamic  pressure. 

Al60  the  aerodynamic  forces  are  applied  to  the  entire  projected 
contact  ellipsoid  area  that  has  penetrated  the  wind  plane  even  if  the 
ellipsoid  is  partially  or  fully  blocked  by  another  ellipsoid.  This 
causes  a  disproportionate  amount  of  force  to  be  applied  in  many  cases. 
This  is  especially  significant  for  the  tor60  segments  where  the 
ellipsoids  sub6tant ival ly  overlap.  More  than  30  percent  cf  three  torso 
segment's  combined  area  is  within  another  ellipsoid,  resulting  in  the 
aerodynamic  force  on  these  segments  being  much  too  large. 

The  original  aerodynamic  forceB  are  applied  to  any  segment  by 
specifying  an  aerodynamic  pressure,  a  boundary  plane,  and  a  contact 
ellipsoid  associated  with  the  segment.  When  the  ellipsoid  penetrates 
the  boundary  plane,  the  wetted  area  is  estimated  and  a  pressure  from  the 
tabular  data,  defining  the  time  dependent  aerodynamic  pressure,  is  used 
to  calculate  the  force  and  torque  that  is  applied  to  the  segment. 

Three  changes  have  been  made  to  the  routines  to  allow  more 
flexibility  in  applying  aerodynamic  forces. 

1.  The  aerodynamic  pressure  can  be  a  function  of  a 
segment'6  velocity. 

2.  A  time  dependent  drag  coefficient  can  be  included  in 
calculating  the  wind  force. 
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3.  An  additional  method  of  calculating  the  wetted  area, 

that  allows  segments  to  be  defined  which  block  the  ellipsoid 
from  the  wind,  has  been  added  as  an  option. 

2.1.  VELOCITY  DEPENDENT  PRESSURE 

To  allow  for  a  velocity  dependent  aerodynamic  pressure,  subroutine 
KINPUT  is  altered  to  read  in  E.6  cards  that  contain  the  specific  heat 
ratio,  the  speed  of  6ound,  and  the  absolute  pressure  for  the  altitude 
which  the  simulation  is  to  represent  along  with  the  definitions  of  two 
segments.  The  aerodynamic  pressure  will  depend  on  the  velocity  of  the 
first  segment  with  respect  to  the  second  segment. 

The  aerodynamic  pressure,  FT,  used  in  WINDY  to  determine  the  aerodynamic 
forces,  is  calculated  from  the  definition  of  dynamic  pressure: 

Ft  *  (i/2)kPa  (v/ci2 

where  k  is  th<  ratio  of  specific  heats 
c  ) 6  the  speed  of  sound 
Pa  is  the  absolute  pressure 

V  is  the  velocity  of  the  first  segment  with  respect 
to  the  reference  segment 

Note  that  FT  is  a  pressure  and  is  multiplied  by  a  wetted  area  in 
subroutine  WINDY  to  determine  the  wind  force  applied  to  a  segment.  FT 
can  be  defined  as  time  dependent  using  the  same  input  cards  as  before, 
or  as  velocity  dependent  by  specifying  a  specific  heat  ratio.  How  FT  is 
applied  to  a  segment  in  subroutine  WINDY,  has  no  functional  dependence 
on  the  methou  used  to  define  FT. 
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2.2  DRAG  COEFFICIENT 


lime  dependent  drag  coefficient  functions  can  be  defined  as  wind  force 
functions  on  the  E.fe  cards.  They  follow  the  same  format  as  the  time 
dependent  wind  force  functions,  although  the  drag  coefficient  is  a 
scalar  quantity  rather  than  a  vector.  Before  the  aerodynamic  pressure 
is  used  in  WINDY,  it  is  multiplied  by  the  drag  coefficient,  Cp. 

^new  =  ^ 

This  can  be  used  to  simulate  the  effects  of  the  drogue  chute  opening  or 
other  events  that  elfect  the  drag.  If  there  is  no  drag  coefficient 
defined,  the  default  value  is  1.0. 

2.3  BLOCKED  WIND 

To  allow  tor  blocking  ot  the  wind,  a  second  method  of  applying  the  wind 
force  has  been  added  to  subroutine  WINDY.  This  involves  projecting  the 
ellipsoid,  to  which  the  aerodynamic  force  is  being  applied,  as  an 
ellipse  to  define  the  wetted  area.  Then  this  ellipse  is  divided  into 
incremental  areas,  whose  center  points  are  checked  for  penetration  of 
the  wind  plane  and  for  blockage  by  other  segments.  Each  area  that 
passes  these  tests  has  the  wind  force  applied  at  its  center  point.  This 
allows  for  overlapping  and  connected  segments.  Since  this  new  grid 
method  can  increase  run  time  significantly,  the  original  method  can 
still  be  U6ed  tor  any  or  all  of  the  segments  to  which  a  wind  force  is 
being  applied,  without  any  changes  to  previous  input  decks. 

Subroutine  WINDY  contains  the  major  changes  that  incorporate  thi s 
new  method  for  applying  the  aerodynamic  forces.  Much  of  the  analysis 
needed  for  this  method  is  based  on  the  derivations  developed  for  the 
VIEW  program  (Ref.  6).  In  WINDY,  after  checking  if  there  is  any 
penetration  of  the  segment  through  the  wind  plane,  and  getting  the  wind 
pressure  from  the  wind  force  functions,  the  program  chooses  a  method  for 


the  wind  force  calculations  depending  on  the  input.  The  original  method 
uses  a  calculated  area  of  the  ellipsoid,  while  the  new  method  allows  for 
blocking  of  the  wind  by  other  segments,  by  using  a  grid  to  determine  the 
area . 


2.3.1  Project  Ellipsoid 

For  the  grid  method,  the  first  step  is  to  set  up  a  coordinate 
system  associated  with  the  wind.  This  viewpoint  coordinate  system  is 
located  at  an  assumed  origin  of  the  wind  with  it's  z-axi6  directed 
towards  the  origin  of  the  inertial  coordinate  system. 

Def ine , 

FT  a6  the  wind  force  vector  (inertial  coord.) 


VP  as  the  origin  of  the  viewpoint  coordinate  system 
(inertial  coord.),  which  is  set  equal  to  -10000FT. 


DVP  as  the  direction  cosine  matrix  tor  transformation  of 
vector  components  from  the  inertial  to  the  viewpoint  coordinate 
sy  stem . 


The  DVP  transformation  is  chosen  such  that  the  X-axis  of  the  viewpoint 
coordinate  system  is  parallel  to  the  X-Y  plane  of  the  inertial 
coordinate  system.  DVP  can  be  calculated  as  follows: 


Let 

ft  =  FT/  i FT l  ,  wh  ich  is  the  unit  wind  force  vector 

f  orm 

XNORM  c  »/f t i ^  +  ft2^,  which  is  the  projected  length  of  ft 
on  the  X-Y  plane. 
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Let 


then  it  can  be  shown  that 

A  ^  ^ 

Xyp  =  <. f 1 2 1  -  f t^J )/XNORM,  i.6  a  unit  vector  normal  to 

A 

ZVp  and  parallel  to  the  X-Y  plane  of  the  inertial 
coordinate  system. 

The  third  unit  vector  can  then  be  obtained  using  the  vector  product 


A  A  A 

Yvp  =  Zvp  x  *vp 

The  transformation  matrix  is  then 

formed  by  placing 

these  unit 

vectors 

in  row  form 

DVP  = 

f  t2/XNORM 
f tLf t3/XNORM 

-f  ti/XMORM 
f t2f t3/XNORM 

0 

-XNORM 

_  ftl 

ft2 

f  t3 

The  contact  ellipsoid  is  projected  onto  a  plane  parallel  to  the  X-Y 
plane  of  the  viewpoint  coordinate  system.  Since  the  viewpoint  is  far 
away  from  the  ellipsoid  and  the  Z-axis  is  nearly  directed  at  the 
ellipsoid,  the  projection  is  assumed  to  be  elliptical.  To  solve  for  an 
ellinse  matrix,  three  radial  vectors  of  the  ellipsoid,  pointing  to  a 
surface  point  that  forms  the  contour  of  the  projected  shadow,  must  be 
determined.  BjD(7-15,M)  is  the  matrix  that  defines  the  surface  points  of 
ellipsoid  M  with  respect  to  its  principal  axes.  First,  this  ellipsoid 
matrix,  BD(7-15,M),  is  transformed  to  the  viewpoint  coordinate  system 
and  designated  as  AM(3,3).  This  is  accomplished  by: 

AM  =  DVP  Dt  BD  D  DVPt 

where  D  is  the  direction  cosine  matrix  that  transforms 

from  t lie  inertial  to  the  ellipsoid  principal  coordinate  system. 
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In  order  to  define  the  projected  ellipse,  three  vectors  are  chosen  which 
lie  respectively  in  the  X-Z  plane,  the  Y-Z  plane,  and  the  (X=Y)-Z  plane 
of  a  coordinate  system  parallel  to  the  viewpoint  system  but  with  its 
origin  at  the  ellipsoid  center.  These  vectors  are  shown  in  Fig.  1,  and 
have  components 


with  R3x  =  &3Y 

As  seen  in  Fig-  1  the  associated  vector  Pj  from  the  viewpoint  to  the  tip 
of  Rj  is  normal  to  the  normal  vector  nj  for  the  point  defined  by  R^  on 
the  ellipsoid.  Therefore, 

nj  •  Pj  =  0  and  iij  =  AM  Rj  1 )  &  2) 

Combining  equations  1  and  2, 

AM  R1  •  Pj  -  Px  •  AM  Rl  ■=  0  3) 

Also  from  the  figure, 

Pj  =  SM  +  Rj  A) 

Substituting  in  eq .  3  for  Pj  from  eq.  4, 

(SM  +  Rt )T  AM  R1  =  0  =  SMT  AM  Rj  +  RTj  AM  Rj  5) 

AM  Rj  =  1  from  the  definition  of  an  ellipsoid.  6) 

then, 

SMT  AM  Rj  =  -1  7) 
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Viewpoint  Y 

vp 


Figure  1  Three  Radial  Vectors 


Subroutine  SOLVE  solves  this  equation  for  the  components  of  R3 .  The 
same  procedure  is  used  to  solve  for  the  other  components  of  52  and  R3. 


To  solve  for  R3,  R2.  and  £3,  we  will  respectively  treat  these  as  Case 
No.  1,  Case  No.  2,  and  Case  No.  3. 


For  all  three  cases  expanding  eq.  7 


(SMX,  SMy,  SMZ) 


AMll 

AMi2 

am13 

AM21 

AM22 

A«23 

AM3i 

am32 

AM33 

R  -  -1 


or 


(SMxAMii  +  SMYAK21  +  S^AM31,  SMxAM3  2  +  SMyAM22  +  SMZAM32, 
SMXAM13  +  SHyAM23  ♦  S>fcAM33)  R  »  -1 


This  can  be  further  reduced  to 


(A1SMX  +  A2SMy  +  A3SMz;)  RXory  +  (A4SMX  +  A5SMy  +  A6SMZ)  Rz  -  -1 


where 


Case  No.  1 

Case  No.  2 

Case  No.  3 

A1  -  AMn 

A1  «  AM12 

A1  *  AM33  +  AMj  2 

A2  *  AN 2^ 

A2  =  AM22 

A2  “  AM2\  +  AM22 

A3  -  AM31 

A3  -  AN32 

A3  «  AM31  +  AM32 
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A4  =  AM13 


A5  *=  AM23 


for  all  cases 


A6  =  AM3 3 


Making  the  farther  substitutions; 


B  A1SMX  +  A2SMy  +  A3Sl^ 
D  *  A4SMX  +  A5SMy  +  A6SMZ 


we  get 


B  RXorY  +  D  Rz  =  -1. 
Solving  for  RXorY 


RXorY  =  '(D/B)RZ  -(1/B). 

Now  the  R  vectors  can  be  written 


-(D/B)RZ  -1/B 
0 

RZ 


R2  - 


0 

-(D/B )RZ  -1/B 


L 


RZ 


r3  ‘ 


(D/B)RZ  -1/B 
( D/B )RZ  -1/B 


Substituting  R  into  eq •  6  and  expanding, 


A7  RXorY2  +  2  A8  RXorY  RZ  +  A6  RZ2  '  1  (9) 

where 


Case  No.  1 

Case  No.  2 

Case  No.  3 

A6  =  AM33 

A6  =  AM33 

A6  =  AM33 

A7  =  AMU 

A7  =  A»22 

A 7  =  AMj  2.  +  2AM]^2  +  AM22 

A8  =  AM13 

A8  =  AM23 

A8  =  AM13  +  AM23 

Substituting  eq . 

8  into  eq .  9, 

A7[-(D/B)RZ  -  1/B ] 2  +  2A8RZ[-(D/B)RZ-1/B J  +  A6RZ2  =  1 
Expanding  and  combining  like  terns, 

IA7(D/B)2+A6-2A8(D/BJRZ2  +  [2A7(D/B2)-2A8(1/B) ]RZ  +  A7/B2-l  =  0 
Therefore, 

RZ  =  -T2  +  Vl22  -  4T1T3  and  RxorY  c  -CD/B)RZ  -  1/B 

2T1 


where 

T1  =  A7(D/B)2  +  A6  -2A8(D/B) 
T2  -  2A7(D/B2)  -  2A8Q/B) 

T3  -  A7(l/B2)  -  1. 
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WINDY  calls  subroutine  SOLVR  which  requires  as  input  variables  the 
values  of  A1 ,  A2,  A3,  A4,  A5,  A6,  A7,  A8,  and  SM  for  each  case. 
Subroutine  SOLVR  returns  corresponding  values  for  ^XorY  an<*  RZ*  The 
three  R  vectors  obtained  satisfy  the  three  dimensional  ellipsoid,  and 
lie  in  the  appropriate  planes. 

These  R  vectors  are  then  projected  as  if  the  viewpoint  was  an 
infinite  distance  away.  Therefore,  the  X  and  Y  components  of  R  are  the 
two-dimensional  projected  vectors,  R2. 


R2^  —  R^  R2y  =  &y 

The  equation  for  the  ellipse  is, 


R2^  AS  R2  =  1  where  AS  - 


ASn 

AS12 

AS21 

Since  there  ate  three  R2  vectors  and  three  independent  components  of  AS, 
AS  can  be  obtained  by  solving  three  equations  simultaneously  which  is 
done  in  subroutine  SOLVA. 


2.3.2  _  Set-up  hnd  Pattern 

To  set  up  the  grid  pattern  for  the  ellipse,  the  majci  and  minor 
axis  vectors  are  needed.  These  vectors  are  found  by  solving  for  the 
eigenvalues  of  the  ellipse  matrix,  AS,  by  imposing  the  condition 

AS  R  =  A  R. 

This  condition  is  true  only  for  vectors  that  represent  the  major  and 
minor  axes  of  the  ellipse. 


"ASi , 

ASl  2 

RX 

A  Rg 

CN 

JO 

< 

- 1 

AS^2 

- 

Ry 

\  Ky 
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Or, 


AS\ >,  ASj  2 


AS | 2  AS22"  \ 


0 


(10) 


The ref ort  , 


=  ASjj  +  AS22  --  \/( AS ^  +  AS22^  “  4  (ASj  j  AS22  ~A^ 2^ )  ^  ^ 

2 

W11I1  the  t-igt'n  vectors, 


1  J 

as12 

KX 

/“T 

1  L  ~AS11 

ry 

111) 


1 

•  2  ~as22 

Rx 

_ 

V  •  2 

1 

rs 

CO 

< 

_Ky. 

These  are  the  major  and  minor  axes  vectors  of  the  ellipse. 


2_.  3  _.  3_Three_  Checks 

With  the  major  and  minor  axes  of  the  ellipse  found,  a  grid  is  laid 
over  the  projected  ellipse  (Fig.  2)  and  each  corner  point  of  the  grid  is 

checked  to  see  if  it  is  in  the  ellipse,  through  the  wind  plane  or  not 

behind  a  blocking  segment.  it  all  are  true,  then  a  wind  force  is 

applied  to  the  incremental  area  (AREA)  shown  shaded  in  Fig.  2. 
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Figure  2  Grid  Overlay 


2.3.3.),  Check  if  Copper  Poinfr  Is  Within  Projected  Ellipse.  $K  is  the 
two  dimensional  vector  to  a  corner  point.  If  &S  Si?<4,  then  551  is 
within  che  ellipse. 
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For  the  next  two  checks,  the  three  dimensional  location  oi  RM  on 
the  ellipsoid  is  needed.  See  Fig.  3.  The  Z  component  of  the  three 
dimensional  vector,  EM,  is  found  by  solving  the  ellipsoid  equation, 


AM^  AMi2  AM33I 

'rmx" 

(RMX,  RMy,  RMZ) 

AMi  2  AM22  AK23 

RMy 

_  AMl3  ^23  ^33 

RMZ 

- 

Expanding  and  solving  for  RMZ 


RMZ  -  -1M2  t  v/lM22  -4TM1  T% 

2TMj  ' 


where 


TM^  ■  AM3  3 

TM2  -  2  (RMX  AM13  +  RMy  AM23) 

TM3  -  RMX2  AMU  +  RMy2  AM22  +  2RMX  RMy  AM^  2  -1 

and  SN1  is  the  vector  from  the  viewpoint  to  RM  on  the  ellipsoid's 
surface  and  is  given  by 

M  -  S>?  +  EM 

with  all  vectors  expressed  in  the  viewpoint  coordinate  system. 


2 . 3 . 3^,2.  Check  if  RM  Is  Penetrating  Wind  Plane.  From  Fig.  4  define 

A 

PL  -  normal  unit  vector  to  wind  plane  (in  the  segment  coordinates 
to  which  plane  is  attached) 
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'  I  1  tpsuiil/P  i  .ini-  I’.iH'l  r.it  inn 
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2. 3. 3. 3  Check  Chat  RM  Ib  Not  Blocked  -v  Any  0£  The  Blocking,  Segments. 

First,  the  ellipsoid  matrices  for  the  segments  that  may  be  blocking  the 
wind  must  be  transformed  to  the  viewpoint  coordinate  system.  This  is 
done  using 

AIi  =  DVPj  £iT  BDj_  D,  DVPj 

Also  Sli  is  defined  as  the  location  of  the  i-th  blocking  segment  in  the 
viewpoint  coordinate  system  (Fig.  5). 


Figure  3  Check  For  Blocked  Wind 

A  line  of  sight  is  defined  as  a  lit-r  normal  to  the  viewpoint  X-Y 
plane  through  the  tip  of  RM,  and  is  used  to  determine  it  the  wind  hits 
the  blocking  segment  before  reaching  RM  (Fig.  5).  The  two-dimensional 
vector,  (X,  Y),  from  the  center  of  the  blocking  ellipsoid  to  the  line  of 
sight  is  used  in  defining  the  point  where  the  Kne  of  sight  enters  the 
blocking  ellipsoid. 
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X  =  SN 1 x  -  S1X 


! 


) 


Y  =  SN 1 y  -  Sly 


The  Z  component 

as  RM2  earlier, 

of  the 

entry 

point,  Bl, 

is  calculated  in  the 

same  manner 

J 

"  Ain 

AIi2 

ai13 

X  ~ 

(X,  Y, 

Bl) 

All  2 

ai22 

ai23 

Y 

=  1 

(12) 

All  3 

ai32 

AI33 

Bl 

It  the  lint-  of  sight  does  not  pass  through  the  blocking  segment,  B1  in 
equation  12  is  a  complex  number.  The  entry  point  could  be  beyond  RM, 
therefore,  the  distances  from  the  X-Y  plane  of  the  viewpoint  coordinate 
system  are  compared.  If  SN12  •  SI2  +  Bl,  then  RM  is  not  blocked. 

Each  of  the  possible  blocking  segments  are  checked  using  this 
method.  If  each  of  tnese  checks  are  true,  then  a  wind  force  is  applied 
to  the  incremental  area  at  RM.  Each  corner  point  is  handled  the  same 
way  and  the  forces  are  totaled  and  added  to  the  III  and  132  arrays. 

2.4.  CHANGES  TO  THE  PROGRAM 

A  new  H  card  is  now  needed  for  the  wind  forces  to  be  output  as 
tabular  time  history®  The  wind  force  applied  to  any  segment  can  be 
output  to  the  tabular  time  histories  in  any  reference  system. 

This  addition  requires  changing  common  block  RSAVE  to: 

COMMON/RSAVE/  XSG(3,20,3),  DPMI(3, 3 , 30 ) ,  LPMI(30),  NSG(9), 

MSG( 20, 9) ,  MCG ,  MCGIN(23,3),  KREF(9) 

The  size  of  NSG,  MSG,  and  KREF  are  increased  to  allow  for  the 
additional  6et  of  tabular  time  histories. 
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Besides  RSAVE,  the  other  changes  require  common  block  WINDER  to  be 
modified  to: 

COMMON/WLNDFR/WTIME(.30),QFU(3, 5),QFV(3,  5),WF(3, 30),  IWIND(30), 
MWSEG(7,30),NFVSEG(b),NFVNT(5),MOWSEG(30,30) 

The  size  of  tWSEG  is  increased  to  include  the  drag  coefficient 
function  number  and  the  number  of  possible  blocking  segments. 

New  variables  that  are  added  to  the  COMMON  block  are: 

WF(.3,  30)  Wind  force  vectors  applied  to  segments 

(.in  local  reference), 

M0WSEGC3G,  30) 

M0WSEG(2I-1,  J)  Segment  identification  number  of  I-.h 
segment  that  can  block  segment  J, 

MCWSEGC2I,  J)  Contact  ellipsoid  associated  with  the 

MOWSEG  C  21—1 , J )  segment. 

In  addition  to  the  common  block  changes,  coding  changes  affected  a 
number  of  subroutines.  Subroutines  WINDY  and  K1NPUT  contain  the 
majority  of  these  changes  and  SOLVA  and  SOLVR  are  new  subroutines.  In 
Volume  3  of  this  report,  the  listing  of  the  ATB-IV  code  has  the  labels, 
WINDOP  or  WINDROT,  in  column  73,  of  all  the  new  or  changed  lines  needed 
for  these  wind  force  options. 

To  use  the  velocity  dependent  pressure,  drag  coefficient,  or  blocked 
wind  the  input  deck  has  to  be  modified.  The  input  description  for 
ATB-IV  16  in  Volume  2  of  this  report  and  describes  these  modifications. 
Note:  Previous  ATB  or  CVS  input  deck6  require  a  blank  card  to  be 

inserted  for  H.8. 
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3.0  ANGULAR  DRIFT  CORRECTION 


The  locked  axes  of  the  joints  in  the  ATI  Model  often  drift  from  their 
original  position  because  of  inherent  inaccuracies  due  to  the  numerical 
integration  process.  The  CHAIN  subroutine  was  written  to  correct  for 
these  errors  after  each  integration  step,  but  the  drift  of  the  locked 
axes  6till  occurred,  especially  during  long  simulations.  The  code 
modifications  described  here  correct  this  drift  and  the  sudden  shifts  in 
the  joint  azimuth  angles. 

3.1  TECHNICAL  DISCUSSION 

The  ATB  Model  has  four  types  of  joints,  they  are: 


1 ) 

Ball  and  Socket  Joint 

2) 

Pin  or  Hinge 

Joint, 

3) 

Euler  Joint, 

and 

4) 

Null  Joint. 

The  Ball  and  Socket  Joint  and  the  Pin  Joint  may  be  locked.  The 
Model  will  unlock  these  joints  when  a  specified  torque  is  exceeded.  The 
Euler  Joint  has  three  axes  which  may  be  locked  or  unlocked  independently 
thus  providing  eight  states  for  the  joint.  The  Null  Joint  is  used  to 
provide  the  option  of  disjoining  sets  of  segments. 

In  the  ATB  Model  constraints  are  imposed  on  the  joints  by  computing 


a  constraining 

torque. 

The  basic  equations 

are : 

II  *1 

+  J2i  Iq 

=  other  torques 

(13a) 

1 2  *2 

~  £2  iq 

=  other  torques 

(13b) 

where  .1^  , 

i2  are 

inertia  matrices  of 

the  adjoining 

segments, 


W|  ,  W2  are  the  angular  accelerations, 
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fil  >  fi.2  are  direction  cosine  matrices, 


P  is  a  projection  matrix  which  depends  on 

the  type  of  constraint,  and 

q  is  the  constraint  torque. 

The  constraint  equation  (which  is  needed  to  solve  for  q)  can  be  best 
derived  by  considering  the  case  of  a  Pin  Joint.  In  this  joint  the 
constraint  is  that  the  pin  vector  in  one  segment  must  coincide  with  the 
pin  vector  in  the  adjoining  segment.  The  pin  vector  defines  the  free 
axis  of  the  join*-  (.the  Pin  Joint  has  only  one  free  axis).  The  equation 
is : 

D^q  =  2^2  h2=h  (14) 

where  ,  h2  are  the  pin  vectors  (1x3  matrices) 

in  the  respective  segments.  These 
vectors  are  constant  in  the 
segments. 

h  Is  the  instantaneous  pin  vector  in 

inertial  reference, 

and 

D^q  >  2^ 2  are  the  transposes  (inverses)  of  the 

direction  cosine  matrices. 

Differentiating  equation  (14)  yields  an  equation  in  the  velocities  w q 
and  V2  thus: 

£T1  («!  x  h^)  =  2^2  x  ^2^  ^5) 
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Differentiating  equation  (15)  yields  an  equation  in  the  acceleration, 
thus : 

DT1  ( W]_  x  hi)  +  (*i  x  (w^  x  h[)) 

(16) 

*  £T2  («2  x  h2)  +  (®2  x  (®2  x  f»2 > ) 

Equation  (16)  is  the  desired  constraint  equation  for  the  accelerations, 
however  it  is  only  of  rank  2  and  we  need  an  equation  of  rank  3  to  solve 
for  the  torque.  This  is  obtained  by  observing  that  the  constraint 
torque  can  have  no  component  on  the  pin  axis,  i.e. 

hTq  =  0  (17) 

Equations  (16)  and  (17)  can  be  combined  into  a  single  matrix  equation  of 
tank  3  by  crossing  equation  (16)  with  h  and  adding  the  term  hhTq.  The 
resulting  equation  can  be  put  in  the  form: 


I  £T1  n  -  £  Et2  w2  +  (i  -  £)q 

=  hT2  ^2  £T2  (i>2  x  w2)  -  hT^  w^  (h^  x  w^ ) 


(18) 


The  projection  matrix  £  is  given  by 

P  =  I,  -  h  hc  where  I.  i  s  the  identity  matrix. 


Equations  (13)  and  equation  (18)  are  the  basic  equationsused  in  the 
ATB  Model  to  form  the  system  equations  which  are  solved  for  the 
constraint  forces  and  torques  and  for  the  linear  and  angular 
accelerations. 


Details  on  the  form  of  equation  (18)  for  a  locked  joint  and  for  an 
Euler  Joint  are  given  in  Volume  1  of  Reference  (4).  The  only 
differences  are  the  form  of  the  projection  matrix  P  and  the  right  hand 
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side  of  equation  (18).  In  particular  for  a  locked  joint  P  is  the 
identity  matrix  and  the  right  hand  side  of  equation  (18)  is  the  zero 
vector  (1x3  matrix). 


3.2  CORRECTION  OF  THE  DRIFT  PROBLEM 

The  p’-ob'  ec.  vith  d»ift  arises  since  the  constraints  are  imposed  on 
the  acceleration  (equation  (18))  and  the  values  of  angular  velocity 
(w's)  and  angular  position  (£’s)  are  obtained  by  numerical  integration. 
Errors  can  arise  because  of  errors  in  the  solution  of  the  system 
equations  and  errors  due  to  the  numerical  integration  process.  Errors 
of  this  nature  are  unavoidable  because  of  the  finite  precision 
calculation  on  a  digital  computer.  ThuB  we  may  find  that  equation  (14) 
and  equation  (15)  are  not  satisfied  to  some  desired  degree  of  precision 
at  some  point  in  the  solution  process.  To  correct  this.  Subroutine 
CHAIN  was  modified. 

Consider  the  following  equations: 

Let  u  =  h2  x  (£2  fiTl  h^),  a  vector  (1x3  matrix). 

If  the  vector  u  is  zero  then  S2  and  h}  are  aligned.  If  u  is  not 
zero  it  is  perpendicular  to  the  pin  vectors  h^  and  h2  and  has  a 
magnitude  which  is  the  sine  of  the  angle  between  the  pin  vectors.  This 
vector  is  used  to  define  the  rotation  operator  that  is  applied  to 
segment  2' 6  direction  cosine  matrix  aligning  the  and  h2  vectors; 

J?2*  =  l  ci  +  uu*7(l  +  c)  -  (u  x)lfi2  (19) 

where  c  is  the  cosine  of  the  angle  between  the  pin  vectors, 

and  (u  x)  is  the  matrix  analogous  to  a  vector  cross  product 

operation. 
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The  velocity  vector  is  also  modified; 


w2*  =  32*  £^1  W1  +  ^2  (b^2  ^2  -  h^i  wi)  (20) 

Subroutine  CHALN  (at  the  option  of  the  user)  modifies  the  direction 
cosine  matrices  JJ  and  the  angular  velocities  w  as  specified  in  equations 
(19)  and  (20).  This  insures  that  the  angular  position  constraint 
(equation  (14))  and  the  angular  velocity  constraint  (equation  15)  are 
satisfied  for  pin  joints. 

3.3  CORRECTIONS  TO  THE  PROBLEM 

The  ATB  model  was  studied  in  detail  to  determine  why  the  above 
procedure  was  not  functioning  properly.  Two  errors  were  found,  they 

are  : 

1.  The  right  hand  side  of  equation  (12),  which  is 
computed  by  Subroutine  SETUP1,  was  being  computed 
before  the  direction  cosine  matrices  were  modified. 

2.  Incorrect  U  vectors  were  used  in  Subroutine  CHAIN 
for  an  Euler  Joint  in  states  4,  5,  or  6. 

Error  1  was  corrected  by  calling  Subroutine  CHAIN  before  calling 
Subroutine  SETUP1  in  Subroutine  DAUX,  and  error  2  was  corrected  by 
correctly  defining  the  h  vectors.  Also  the  code  to  correct  for  drift 
was  removed  from  subroutine  CHAIN  and  put  into  a  new  subroutine  called 
DRIFT.  This  required  that  the  dimensions  of  the  HIR  array  and  the  CONST 
array,  which  are  in  COMMON/CEULER/,  be  changed.  The  new  arrays  are 
HIR(3,3,90)  and  CONST(5,30).  This  change  was  made  in  all  the  routines 
that  included  this  COMMON. 

Subroutines  EJOINT,  INITIAL,  and  UPDATE  were  modified  to  store  the 
variables  needed  for  the  new  drift  routine. 
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Listings  of  the  new  DRIFT  subroutines  and  of  the  changed 
subroutines  are  in  Volume  3  of  thi6  report.  New  or  changed  lines  are 
labeled  with  JDRIFT  starting  in  column  73. 
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4.0  EDGE  EFFECT  OPTION 


In  the  past,  problems  often  developed  in  ATB  simulations  when  an 
ellipsoid  came  in  contact  with  a  plane  near  the  plane's  edge.  If  part 
of  the  ellipsoid  contacted  the  plane  edge  but  the  center,  of  the 
crosB-sectional  ellipse,  containing  the  area  cut  by  the  plane,  did  not 
lie  within  the  plane  boundaries,  then  no  force  was  applied  as  if  no 
contact  had  occurred  at  all.  However,  as  60on  as  this  center  moved 
within  the  boundaries  of  the  plane,  a  full  contact  force  was  applied. 

The  planes  had  to  be  adjusted  and  modified  frequently  to  avoid 
instantaneous  jumps  in  force  when  contact  occurred  at  object  corners. 
These  new  routines  have  been  developed  to  solve  this  edge  effect 
problem.  In  particular,  use  of  the  new  edge  effect  option  insures  that 
a  contact  of  a  plane  with  an  ellipsoid  will  not  be  ignored  and  that  a 
smaller  force  will  be  applied  when  only  part  of  the  contact  area  is 
within  the  plane  boundaries.  Also  an  option  was  added  allowing  a  force 
to  be  applied  when  the  ellipsoid  has  completely  penetrated  the  plane,  if 
the  edge  effect  option  is  not  used.  Another  alternative  for  improving 
contact  force  calculations  is  to  use  a  hy perel 1 ip6cid  to  describe  the 
surface.  This  option  is  described  in  Section  Seven. 

4.1  NEW  SUBROUTINES 

4 . L . I  Subroutine  PLELP 

Subroutine  PLELP  computes  the  point  of  maximum  penetration  of  an 
ellipsoid  associated  with  segment  m  intersecting  a  plane  associated  with 
segment  n.  Previously  the  point  of  maximum  penetration  was  projected 
onto  the  plane.  If  this  projection  tell  outside  of  the  boundaries  of 
the  plane,  the  contact  was  ignored.  A  five  way  option  has  been  added  to 
the  routine.  The  choice  is  made  by  the  user  by  inputting  an  additional 
integer  on  the  F.l.B  -  F.l.N  cards.  This  integer  is  stored  in  the 
twenty  third  location  of  the  TAB  array  associated  with  the  contact.  The 
options  are: 
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TAB  ( NT-t-22  )  , 


0,  call  new  edge  effect  routine  PLEDG,  no 

force  is  computed  for  complete  penetration, 

=  0,  use  standard  finite  plane  test,  no  force  is 

computed  for  complete  penetration, 

=  -1,  treat  plane  as  infinite  (bypass 

edge  test),  no  force  is  computed  tor  complete 
penetrati on, 

=  -2,  treat  plane  as  infinite  (bypass  edge  test),  a 
force  is  computed  for  complete  penetration, 

'  -2,  use  standard  finite  plane  test,  a  force  is 
computed  for  complete  penetration. 


Equations  used  PLELP. 


Let:  (Fig.  6 


Fl 

A 

T 

E 


location  of  the  reference  point  of  segment  m, 

( inertial  system) 

location  ot  the  reference  point  or  segment  n, 
(inertial  system) 

offset  of  tne  ellipsoid,  (inertial  system) 
first  reference  point  for  the  plane,  (inertial 
sy  8 1  em  ) 

unit  exterior  normal  of  the  plane, 
ellipsoid  matrix. 


Then  the  equations  are: 


*nc  “  zm  +  -  Zn  -  PI  ,  vector  from  PI  to  center  of 

el  1 ipsoid . 
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bet  =  -  l^Xjjc  ,  distance  from  plane  tc  center  ot 
ellipsoid  '.positive  if  center  or 
ellipsoid  hat;  penetrated  the  plane.'. 

bte  =  -sqrt(T^  E-^  T),  the  component,  normal  to  the  plane, 

of  the  vector  lion  the  point  of 
maximum  penetration  to  the  ellipsoid 
ante r  . 

—  A 

Rji  =  £  1  T  /  bte  ,  vector  from  center  of  ellipsoid  to 

point  ot  maximum  penetration.  If 
the  eng-.-  rool.oe  is  used  this  vector 
is  computed  by  PLEDG  as  the  location 
o.  the  ceatrcu)  ct  the  common  area  of 
into  r set 1 1  on . 


p  -  bet  -  bt«.  ,  penetration.  It  the  edge  routine  is 
i:m  i,  t  tu  penetration  is  computed  by 
i'LhiJa  hr-  tie  point,  on  the  ellipsoid 
'e.ic  along  t  i"  :..rn.ai  t.o  the  plane) 
tie  .  e >-.tr. 

amr  •  1  hie  *  pert-Pct  vt  s  tc  determine 

it  silt;  t ; .  .as  occurred 

■  frt  ■  .’lit  . 1  i  t  ;  i'i  •  IX> ) . 

Note  that  it  amr  is  nega .  •  a,.-  1  u  sn  -1  doesn't  intersect  the 

plane.  If  amr  is  zero  the  ihoas'  ;s  lament  to  the  plane.  The 
current  logic  in  Subroutine  i'L  id. !  ignores  'he  uitrii  act  unless  amr  is 
greater  than  zero.  Thi  s  has  .  .*  eite.t  el  Mopping  the  contact  once  the 
ellipsoid  has  fully  penet rat -tl  the  piano. 

Rl^  ~  R^  i  (i^  ,  location  ot  the  poira  ut  maximum 
penetration  relative  to  the 
reference  point  of  segment  m. 
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Rln  =  **m  +  +  _  Za,  location  of  the  point  of 

maximum  penetration 
relative  to  the  reference 
point  of  segment  n. 

4.1.2  Subroutine  PLEDG 

If  the  ellipsoid  intersects  the  plane,  the  figure  of  intersection 
will  be  an  ellipse.  Subroutine  PLEDG  is  used  to  determine  if  this 
ellipse  has  any  common  area  with  the  finite  plane. 

Z  ■=  aU  +  bV  is  a  point  in  the  plane  where 

U,  V  are  the  vectors  defining  the  boundary  of  the 
finite  plane,  and  (U  =  P2  -  PI,  V  '  P3  -  PI, 
where  PI,  P2,  P3  are  the  vectors  defining  the  plane), 
and 

a,  b  are  scalars.  Let 

W  =  (bet/btelR^,  be  the  vector  from  the  center  of  the 
ellipsoid  to  the  center  of  the  cross - 
s  ctional  ellipse. 

d^,  d£  is  location  of  the  center  of  the  ellipse  from 
the  reference  point  PI  in  U,  V  coordinates. 

Then,  the  equation  of  the  ellipse  is 

(Z  -  W)T  E  (Z  -  W)  -  1 . 


This  may  be  written  as 


a^ell  +  2abej2  +  b^e22  “  p  2,  white  (21) 

eu  =  UTflU,  e12  *  UTE^,  e22  =  VT&V,  and  p  2  ---  1  -  WT0? 
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The  extreme  values  of  s  Etui  b  for  t:  r»  i-  ellipse  are: 

(Fig.  1) 

aa2  =  sqrt  ft'22  •  2/  deft,1  >  largest  value  or  a, 

ba2  =  “  eI2aa2/e22>  val «r  3  at  a32, 

aa^  ~  -  aa2  i  smallest  value  of  a, 

bai  -ba2  >  cf  b  at  , 

bb2  =  sqrt  (ep^.  2/dclt/  ,  largest  value  of  t, 

ab2  =  _e12Db2''ell>  valee  ot  u  a!  r  o  2  • 

bbi  -  "bb2,  smallest  Value  ol  o, 

aoL  _3d2  ’  value  oi  a  ar.  b:.- , 

where 


Computation  of  the  intersection  of  the  ellipse  and  the  finite  plane  (a 
parallelogram)  i6  performed  as  follows: 

Def ine 

=  max(aa^,  -dj),  and 
%ax  =  ®in(aa2,  1  -  d1). 


If  is  greater  or  equal  to  a^*,  the  ellipse  has  no  common  area 

with  the  finite  plane.  (It  lies  entirely  to  the  right  of  the  finite 
plane  or  entirely  to  the  left.) 

Def ine 


b„,in  =  max(bbl,  -  do),  and 
^nax  =  mni(bb2,  1  -  d2). 

^  ^nin  is  greater  or  equal  to  h,^*,  the  ellipse  has  no  common  area 
with  the  finite  plane.  (it  lies  entirely  above  the  finite  plane  or 
entirely  below.) 


It  ^nn  18  greater  than  bbj,  the  lower  boundary  of  the  plane 
intersects  the  ellipse.  The  corresponding  values  of  a  may  be  calculated 
from  the  ellipse  equation  21  as 


afmr,  ~  "bni in  e 1 2 / e 1 1  -  8qrt[>  2/ejj-  del t( „/ ej j )2 ] 
af  pn  =  -bmin‘-12/ni  +  sqrtl,  2/e^-  del  t  ( bmln/ eu  ) 1 J 


Similarily  it  bmax  is  less  than  bb2>  the  corresponding  values 

ot  a  are 
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afmx  ^naxe12'eli  aqrti  ,  2/e^A  flfci  t  ( ei  i  ^  i 


afpx  =  -,Doaxe12/ elL  +  b£lr!:i  p  2/eu-delt(bmax/e11  j2] 


h 


l  --in.  -n  \t'i  ,i  bound  . i  i  • 


Peternng  lo  Figure  8,  the  most  general  intersection  ot  the  plane 
and  the  ellipse  consists  ot  three  sections  comprising  an  upper  boundary 
and  three  sections  comprising  a  lower  boundary.  The  computation  of  the 
area  and  the  centroid  is  done  in  subroutine  PLRhA  using  the  abscissa,  a, 
as  the  independent  variable.  The  abscissa  of  the  sections  are 
determined  as  follows: 


Section 


I  left  upper  arc  of  ellipse: 

initially  a^  =  a^  =  amin*  then 

if  fSnax  >  bal>  %2  =  ®in^afmx»  ^ax^ 

if  ^max  11  bal>  an<i  ^min  >  bal» 
ajul  -  nax(a£mn,  ^in^ 

II  upper  straight  line  boundary: 

afp  =  min(a£pX,  a^ax) 

afm  =  max(a£mx,  ajj,£n) 

III  right  upper  arc  ot  ellipse: 

initially  axl  *  ax2  =  amax>  then 

if  ^max  >  ba2>  axl  =  max  ^afpx  >  amin^ 

if  hmax  '  ba2>  anc*  hmin  >  ba2> 
ax2  =  mln  ^afpn>  amax^ 

1*  left  lower  arc  of  ellipse: 

initially  bml  -  b^  =  amin>  then 

if  hffiin  bal»  *^2  =  min  iafmn> 

if  hmin  bal>  and  b^x  '  bai  > 
hml  =  max(a£mx,  a^j^J 
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II'  Lower  straight  line  Doundary  : 

bfp  -  minla^pjj,  a^ax-1 

b£m  =  -«<*£*„-  pm  i  u  ] 

III'  right  lower  ire  of  ellipse: 

initially  bxl  -  bx2  =  a,^,  then 

lt;  ^mm  baz>  bxi  -  max(afpr;,  a^j,) 

*■*  ^min  ba2>  anc*  ^max  ba2> 
bx2  =  min(afpx>  a,,,^) 

With  the  above  logic,  it  the  lower  aoscissa  oi  any  section  16 
greater  than  or  equal  to  the  upper  limit  the  section  does  not  exist.  If 
uc  section  ol  the  upper  boundary  exists  then  there  is  r.c  common  area  and 
the  routine  exits. 

if  it  is  determined  that  a  cutunon  area  exists,  subroutine  PLREA  is 
called  to  compute  the  centroid  of  the  common  area. 

Next,  tne  vector  from  the  center  of  the  ellipsoid  to  the  centroid 
of  the  common  urea  is  computed  as 

Rjj,  =  a(_UbtV  +  Vi,  where 

ac,  bc  are  the  abscissa  and  the  ordinate  of  the 
centroid  and 

W  is  the  vector  from  the  center  of  the 

ellipsoid  to  the  center  of  the  ellipse 
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Finally,  the  penetration  parameter,  p,  is  computed  as  the  point  on 
the  ellipsoid  below  the  centroid.  This  is  done  by  using  the  equation  of 
the  ellipsoid 

( -  p I >  -  1,  where  ^  i6  the  normal  to  the 

plane. 

This  equation  is  a  quadratic  in  the  parameter  p  and  may  be  readily 
solved  with  the  quadratic  formula.  The  penetration  parameter  p  is  the 
positive  solution  of  this  quadratic  equation. 

4.1.3  Subroutine  PLREA 


This  routine  computes  the  common  area  and  its  centroid.  Since  the 
equation  of  the  ellipse  (equation  21)  is  a  quadratic  the  integration  can 
be  done  in  closed  form.  The  abscissa  is  used  as  the  independent 
variable.  For  a  given  value  of  the  abscissa,  a,  the  ordinate,  b,  is 
computed  from  equation  21  a6 

b  -  -aej2/e22  +  sqrt[r/e22  ~  d(a/e 22^J>  where  (22) 

r  =  2  and  d  =  delt,  are  already  computed  by  PLEDG. 

The  area  is  computed  by  adding  the  area  contributed  by  each  ellipse 
or  straight  line  section.  For  the  ellipse  sections  incremental  areas 
are  obtained  by  integrating  equation  22  for  each  arc: 

/  bda  =  h2[t  +  sm(t)  cos(t)  -  f  sin^(  t  )/g]/(  2g) 

For  the  straight  line  portions  the  incremental  area  is  ab. 

The  abscissa  of  the  centroid  is  determined  by  summing  the 
contributions  from  each  section  and  then  dividing  by  the  total  area. 

For  the  ellipse  sections  the  contributions  are  determined  by  integrating 
equation  22  times  a: 

/  bada  =  h^l-cos-^(t)  -  f  sin^(  t  )/g  ]  /  ( 3g^) 


37 


For  the  straight  line  portions  the  contribution  is  ba2/2. 

Similarly  the  ordinate  oi  the  centroid  is  determined  by  summing  the 
contributions  from  each  section  and  divining  by  the  area.  The  ellipse 
contributions  are  equal  to  one  half  the  integral  of  the  square  of 
equation  22: 

1/2  >  b2da  -  h3  [  2t  cos-'(  t  )  *  3gsir.lt)  +  ( t  2  -  g£)  sm3(  t  )/g]/(6g2) , 
and  the  6tr2ighc  line  contributions  are  ab^/2. 

Whert  sid(‘J  -  ag/n, 

f  *=  t12/e 22. 
g  =  sq  r t ( d  ,'/ 1 22>  and 
h  =  sqrtl  r/f;22)  • 

The  routine  combines  these  contributions  to  compute  the  centroid 
and  the  area.  The  true  area  is  the  area  competed  by  thi 6  routine  times 
the  magnitude  of  the  cioss  product  of  the  U  and  V  vectors.  The  true 
area  is  not  computed  since  it  is  not  currently  used. 

4.2  Modification  of  ctk :-;k  routines 

In  the  development  of  the  edge  effect  routine  it  was  convenient  to 
have  the  vectors  U  and  V,  which  are  the  sides  of  the  parallelogram  of 
the  finite  plane,  available.  Hence  the  dimension  of  the  plane  array  in 
COMMON /CNTSRF/  was  changed  from  17  to  2.4.  This  change  was  made  in  all 
routines  that  included  this  COMMON.  The  new  format  of  the  plane  array 
is  given  in  the  table  below. 
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Table  1 


Format  of  PL  Array 


subscript 


Description 


1,  2,  and  3 
4 

5,  6,  and  7 
8,  9,  and  10 

11 

12 

13,  14,  and  15 

16 

17 

18,  19,  and  20 
21,  22,  and  23 
24 


unit  exterior  normal,  T,  to  plane, 

$TPl 

PI 

A 

Up,  unit  vector  perpendicular  to  side 

P3  -  PI 

UpTpl 

UpT(P2  -PI) 

* 

Vp,  unit  vector  perpendicular  to  side 

P2  -  PI 

VpTPl 

VpT(P3  -  PI) 

P2  -  PI,  o’  vector 
P3  -  PI,  7  vector 
not  currently  used 


PI,  P2  and  P3  are  the  points  defining  the  plane. 


Subroutine  ROTATE  was  modified  to  rotate  the  proper  components  of 
the  new  PL  array. 


Subroutine  FINPUT  was  modified  to  allow  the  input  of  the  parameter 
used  to  select  the  edge  effect  option  in  subroutine  PLELP. 

Subroutine  EQU1LB  was  modified  to  use  the  new  PLELP  instead  of 
using  a  shortened  version  of  the  subroutine. 


In  making  changes  to  PLELP,  a  number  of  improvements  were  developed 
to  sections  of  code  that  are  in  both  PLELP  and  SEGSEG.  Therefore, 

SEGSEG  was  changed  to  incorporate  these  improvements. 
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Listings  of  th->  changed  subroutines  and  new  PLEBG  and  PLSEA  are  in 
Volume  3.  New  or  changed  lines  have  EDGE  in  column  73. 

To  use  this  new  option  changes  .0  the  1.1  cards  are  required.  The 
changes  are  described  in  the  input  description  in  Volume  2. 
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5.0  MULTI-AXIS  ANGULAR  VEHICLE  DISPLACEMENT 


ATB-II  and  subsequent  versions  of  the  ATB  Model  allow  the  user  to  input 
a  specified  motion  for  predesignated  segments.  Thus  input  may  be  in  the 
form  of  positions,  velocities  or  accelerations.  However,  position 
information  for  angular  orientation  may  be  specified  only  for  one  axis, 
i.e.  yaw,  pitch  or  roll.  In  order  to  remove  this  restriction  it  was 
necessary  to  develop  a  technique  that  will  allow  general  angular 
orientation  data  to  be  u6ed. 

5.1  MATHEMATICAL  DEVELOPMENT 

Given  the  yaw,  pitch  and  roll  angles  of  a  body  at  specified  points  in 
t  me : 

Let  TM(n)  be  the  time  at  point  n,  where  n  =  1,  N 

(these  time  points  need  not  be  equally  spaced), 

A1 ( n, 1 )  be  the  yaw  angle  at  point  n, 

Al(n,2)  be  the  pitch  angle  at  point  n, 

A1 ( n, 3 )  be  the  roll  angle  at  point  n,  and 
N  number  of  data  sets. 

The  angular  velocity  V  may  be  calculated  from  the  quaternion  product 

V  =  2  q*  q  (23) 

and  the  angular  acceleration  from, 

A  =  2q*  i}D  (24) 

where  V  is  a  column  matrix  (3x1)  with  components  V(i),  l  =  1,3 
A  is  a  column  matrix  (3x1)  with  components  A(i),  i  =  1,3 
q  is  the  quaternion  representing  the  angular  orientation, 
q  is  the  conjugate  of  the  quaternion  q, 
q  is  the  time  derivative  of  the  quaternion  q,  and 
q  is  the  second  time  derivative  of  the  quaternion  q. 
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The  quaternion  q  may  be  represented  as  a  column  matrix  with 
components  q(i),  i  =  1,4.  It  is  convenient  to  consider  the  first 
component  q(L)  as  a  scalar  and  the  other  three  as  a  vector,  i.e.  a  3x1 
matrix  hence  let 

c  =  q  ( 1  ) 

u  =  a  3x1  matrix  with  components  q(i),  i  =  2,4  and  let 

e  =  q( 1 ) 

v  =  a  3x1  matrix  with  components  q(i),  l  =  2,4 

e  =  q  (  i  / 

v  =  a  3x1  matrix  with  components  q(i),  i  “  2,4 


The  quaternion  q  may  be  determined  from  the  direction  cosine  matrix 
representing  the  angular  orientation  of  the  body.  The  direction  cosine 
matrix  is  obtained  from  the  yaw,  pitch  and  roll  angles. 


D  =  (c^  -  u^u)I  +  2  u  u^  -2c  (u  x) 


where 


D  is  the  direction  cosine  matrix  (3x3), 

uT  is  the  transpose  of  the  matrix  u, 

(u  x)  is  the  matrix  representing  the  cross  product  operation 


and, 


I  is  the  identity  matrix. 


From  the  above  equation  one  has 


jT  -  £  *  4  c  (  u  x)  and 

trace(J))  =  3  c^  -  3^0  =  4  c^  -  1  since  c^  +  =*  1. 


Therefore  Q ( 1 ) 
Q(2) 
Q(3) 
Q(4) 


c  =  0.3  SQ ET( D( 1 , 1 )  +  D( 2, 2)  +  D(3,3)  +  1), 
u(l)  =  ( D( 2, 3 )  -  D(3,2))/(4  c), 
u( 2)  •=  (D(3, 1)  -  D(  1 , 3)  >/ (4  c), 
u(  3  )  =  ( D(  1 , 2)  -  D(  2, 1 ) )/  (4  c). 
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It  should  be  noted  that  the  magnitude  of  the  u(i)  may  be  computed  from 
the  formula  u(i)  =  SQRTC(D(i,iJ  +  l  -  c^)/2).  This  may  be  u6ed  in  the 
special  cases  where  c  =  0.  The  signs  of  the  u(i)  may  be  resolved  by 

examining  the  signs  of  the  off  diagonal  terms  of  +  fi. 

Assuming  that  these  angles  are  continuous  in  time,  c  and  u  can  be  fitted 

with  the  Spline  Routine  that  is  already  part  of  the  ATB  Model.  The 
Spline  Subroutine  produces  a  set  of  functions  SP(i,j,k)  for  a  cubic 
spline  which  preserves  the  given  values  of  the  quaternion: 


sp(i, j,k) 

=  TM( j  ) 

,j  =  1,N,  k  = 

1,4 

SP( 2, j , k) 

=  qfj.k) 

,j  =  1,N,  k  = 

1,4 

SP(3, j , k) 

=  ql(j.k) 

,j  -  1,N,  k  = 

1,4 

SP(4, j , k) 

=  q  2  C  j  ,  k ) 

,j  =  1 , N ,  k  = 

1,4 

SP( 5, j  ,  k) 

=  q3(j,k) 

, j  =  1 ,N,  k  = 

1,4 

where  ql,  q2  and  q3  are  the  linear,  quadratic  and  cubic  terms  determined 
by  the  spline  routine.  These  were  determined  such  that  the  first  and 
second  derivatives  of  the  piecewise  cubics  are  continuous  at  the 
specified  time  points  and  the  sum  of  the  squares  of  the  changes  in  the 
third  derivative  at  these  points  is  minimized. 

Using  these  spline  functions  as  interpolating  functions,  values  of  the 
quaternion  terms  may  be  determined  at  intermediate  time  points  by  the 
formula 

BU.lt)  =  SPU.m.k)  +  X  (SP(3,m,k)  +  X  (SP(4,m,k)  + 

X  SP(fi,m,  k)l) 

where  B(t,k)  is  the  interpolated  value  of  quaternion  term  k  at  time  t, 

X  =  t  -  SP(l,m,k)  and  m  is  selected  such  that, 

X  is  positive  or  zero  and  t  -  SP(l,m,k+l)  is  negative. 

The  values  of  the  time  derivatives  of  q  are  estimated  using  the 
derivative  of  the  6pline  interpolating  formula  at  each  time  point, 


q(  k)  =  SP(3,k,3)  +  X  <2.0SP(4,k,U  -r  3.0  X  SP(5,k,L) 


qU)  -  2.0  SP(4,k,lu)  +  0.0  X  SP(5,k,L) 

The  values  of  q,  q,  and  q  are  then  used  in  equations  23  and  24  to 
calculate  the  angular  velocity  and  acceleration  at  the  specified  time 
points  ; 


V  =  2 lev  -  eu  +  vxu  i 


in  detai 1  : 


v(  i ;  -  2  ( 

C 

V  (  1  1 

~  e 

u(l) 

f 

v  (  2) 

u(3;  -  v  k  3 ) 

u(  2 ) 

V \2 j  -  2  ( 

c 

v(2; 

-  e 

U 1  2 , 

v  v  3  J 

U 1  J  -  V  (  1  /' 

u(  3 ) 

V  i  3  !>  =  2  ( 

c 

v  (3 : 

0 

v(3j 

■f 

v  (  1  ) 

j  (  2  ,*  -  v  1  2  .) 

u(l) 

A  -  2i  <:v  - 

4u  +  V 

X  u 

) 

detai  1 

A  ( 1  )  —  2 1 

c 

v  'v  i  ; 

-  e 

nil) 

*- 

v  i  2 ) 

u ( 3  )  -  v(3) 

uU)) 

a;  2  =  2 t 

c 

v  ( 2 ) 

-  V 

u(2J 

t 

v  ( 3  ) 

■i(l)  -  v(l) 

u(  3 ) ) 

a  (  j  .)  =■  ;■(, 

c 

V  '  $  j 

-  e 

u 1  3 , 

4 

v(  n 

u(2l  -  H2) 

u  ( 1 )  ) 

The  angular  acceleration  is  saved  and  used  to  prescribe  the  vehicle 
motion  as  it  is  lor  the  other  cases  when  1  he  angular  velocity  or 
acceleration  is  input. 

3.2  CHANGES  TO  THE  PROGRAM 

The  above  calculations  are  done  in  subroutine,  VLNPUT.  A  new 
subroutine,  QUA!  has  been  added  to  calculate  a  quaternion  from  the  yaw, 
pitch,  and  roll  angles.  Both  VLNPUT  and  QUAT  are  listed  in  Volume  3. 
The  changed  lines  in  VLNPUT  are  labelled  with  JTF984  in  column  73. 
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6.0  SLIP  JOINT 


The  ATB  Model  insures  that  joints  do  not  pull  apart  by  imposing  a 
constraint  force.  A  new  feature  that  has  been  added  to  the  model  is 
termed  a  "slip  joint".  The  slip  joint  allows  one  segment  of  a 
particular  tree  structure  to  move  linearly  along  a  prescribed  axis  with 
respect  to  the  adjoining  segment. 

6.1  EQUATIONS 

First,  we  look  at  the  equauicns  ror  a  two  segment  model. 


The  linear  equations  o*  motion  are: 

ml  *1  *  *  *  ^li 
m  _  x  _■  - U1: 

subject  to  Hit  constraint  t-.uari.on: 
X]  +0]  ^7  t  -  >;  j  ♦  I';  '  T  2 

where 


(25) 


(  26) 


m  j  ,  m  2 

*1>  *2 


are  the  masses  o;  segment  1  and  segment  2, 


are  ttie  locations  of  the  centers  of  gravity  (inertial 
s  y  st  em  ) , 


*1*  *2 


rj,  r2 


are  the  linear  accelerations  (inertial  system), 

are  the  locations  of  the  joint  connecting  the  two 
segments  (local  system), 


U 1 1 ,  U22  are  the  external  t  arces  acting  on  the  segments 
f  inert  'at  system),  arid 

i  is  the  constraint  lout  at  the  j^int  (inertial 

sy  stem  J . 


t 

< 

I 

I 


for  the  standard  joint .  r|  ar.u  7 
system;.  of  the  rtt/.ti'  iv<  segiren 
For  the  slip  joint,  equation  26 


are  fixed  in  the 
0 ,  1 . p .  ! u e  joint 
s  rep  la  it’d  by  the 


local  r ef  erence 
does  not  all  apart, 
equation 


x  1 


^21? 


where  q  ; *  ui t placement  of  the  joint  relative  to  the  fixed  position  rj. 
The  vector  q  1 s  constrained  to  lit  along  a  vector  fixed  in  segment  1 . 


Thus,  for  the  displacement  or  iht  joint 

q  ~  a  h  (28) 

vhere  a  is  a  scaisr  i initialized  to  0)  and 
h  is  a  vi  a.-r  fixed  m  segment  i  • 

v  i:  is  selected  as  the  z  axis  of  the  local  joint 
reference  syst«r,  in  segment  1). 

In  order  to  solve  the  system  equations,  equation  T  '  rust  be  doubly 
differentiated  to  obtain  an  equation  for  the  accelerations.  Performing 
this  differentiation  produces  a  term  involving  the  acceleration  of  the 
scalar  a.  Since  this  is  not  known,  it  must  be  eliminated.  This  is  done 
by  eliminating  the  component  of  the  differentiated  equation  27  that  is 
parallel  to  the  rector  h  resulting  in  a  vector  equation  of  rank  2.  To 
produce  a  complete  set  of  equations  an  equation  which  states  that  the 
constraint  force,  f,  have  no  component  along  the  tree  axis,  h  is  added 
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to  the  sy  stem  equations.  This  process  is  similar  to  the  process  used 
for  joint  constraints  as  described  in  Volume  1  of  reference  4. 


6.2  IMPLEMENTATION 

To  define  a  joint  as  a  slip  joint  the  user  is  asked  to  provide 
additional  information  on  the  first  card  of  the  two  card  set  B.3.A1  - 
B.3J1.  On  this  card  the  user  specifies: 

JOINT(J)  as  before, 

JS( J )  as  before, 

JNT(J)  as  before, 

IPIN ( J  )  as  before  plus  these  new  options; 

=  3  slip  join'  with  complete  angular  freedom 

(same  as  IPIN  =  2  for  angular  motion), 

=■  t>  iip  joint  with  pin  as  y-axis  of  joint 

'■.same  as  LP1N  -  1  for  angular  motion), 

(.program  Will  use  the  flexural  spring 
characteristics), 

-  7  slip  joint  with  pin  as  z-axis  of  the  joint  in 

segment  JNT(J), 

(program  will  use  the  torsional  spring 
charact  er l sties), 

Negative  numbers  may  be  used  to  indicate  that 
the  joint  is  initially  locked,  however,  if  a 
negative  n ud. b*.  r  is  used  ,o  indicate  an 
initially  locked  slip  joint,  ISLIP  cannot  be  0 
(see  1ST  IP  below). 


bR(I,2*J-ll  as  beiore, 
SK( I, 2*J )  as  before, 


ISLIP 


new  parameter, 


=  1  joint  is  a  slip  joint, 

=  0  standard  joint  or  ua Locked  slip  joint, 

=-l  joint  is  a  slip  joint  which  may  be  locked  or 

unlocked  tor  angular  motion  (depending  on 
or  -  IPjN),  but  is  locked  i or  linear  motion. 

Cl  new  paiaueitt,  a  negative  number  for  the  value  of 

unlocking  tore-,  in;  tension, 

C2  new  parameter,  a  positive  n-jxbc-r  lor  the  value  of 

unlocking  force  for  rent;,  rersiou. 

If  :  d’ia  u4  ..re  be  tii  zero  the  joint  v  i  i  1  not  unlock 
f  .  ■  r  B  !  ’  up  a 

The  slip  joints  allow  tot  J  +  I  _jc:r.t  coord. date  system  to  move  along 
tite  z-axis  ot  the  JKT(  J )  ;oint  coord?,  r.a.e  ryslerc.  Tension  is  movement 
along  the  positiv-  z- axis  and  compression  along  the  negative  z-axis. 

When  the  slip  joint  is  tree  to  slip,  spring  and  camper  forces  may  be 
introduced  on  the  slip  axis  by  using  the  spring  camper  option  in  the 
model.  The  spring  damper  option  iJ  speciiii c  on  the  D  8  cards.  The 
coordinates  ot  the  attachment  points  should  be  selected  as  the  location 
of  the  joint  in  the  respective  sfguie  nts. 

6.3  CHANG ES  TO  PROGRAM 

The  dimension  of  parameter  SR  in  COMMON/ DESCRF/  has  been  changed  to 
SU(4,60).  For  joint  J  the  value  of  ’a’  (see  equation  28  above)  is 
stored  in  SR(4,2*J-1)  and  the  value  of  the  time  derivative  of  'a'  i6 
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stored  in  SR  (4,2*J).  These  values  are  initialized  to  zero  by 
Subroutine  BINPUT  and  are  computed  by  Subroutine  CHAIN.  The  change  was 
made  in  all  routines  using  this  COMMON  block. 

A  new  matrix  All  has  been  added  to  COMMON/ CMATRX/ .  The  dimensions  are 
All(3,3,30).  Refer  to  section  4.7  of  reference  4  to  see  how  this  matrix 
is  U6ed  in  the  system  equations  (the  matrix  Bll  is  the  transpose  of 
All).  Matrix  All  is  computed  in  Subroutine  SETUP1.  This  change  to  the 
COMMON  block  was  made  in  all  routines  which  use  this  COMMON  block. 

Subroutine  BINPUT  was  modified  to  input  the  parameters  1SLIP,  Cl,  and  C2 
(see  the  above  section  on  implementation).  The  value  of  1SLIP  is  stored 
in  the  IEULER  array  and  values  of  Cl  and  C2  are  stored  in  the  CONST 
array.  These  arrays  are  in  COMMON/ C EULER/ . 

Subroutine  CHAIN  was  modified  to  compute  the  value  of  'a'  and  its 
derivative  for  a  slip  joint.  The  values  of  SEGLP  and  SEGLV  for  segment 
J+i  (joint  J)  are  adjusted  to  insure  that  the  slip  is  on  the  prescribed 
axis  (axis  h  equation  28  above). 

The  DAUX  subroutines  were  modified  to  introduce  the  matrix  All  into  the 
system  equations  and  account  for  values  of  IPIN  of  3,  6,  or  7. 

Subroutines  DHHPIN  and  DR1ET  were  modified  to  account  for  the  new  values 
of  IPIN. 

Subroutine  PDAUX  was  modified  to  allow  for  the  integration  of  the  linear 
position  and  velocity  of  segment  J+l  if  joint  J  is  a  slip  joint  and  is 
free  to  sLide. 

Subroutine  ROTATE  was  modified  to  accommodate  the  slip  joint. 

Subroutines  RSTART  and  SEARCH  were  updated  to  account  for  the  change  in 
the  dimension  of  SR  and  the  new  matrix  All. 


Subroutine  SETUP!  was  modified  to  compute  the  All  matrix  and  the  VI 
array  for  the  joints. 

Subroutine  SETUP2  was  modified  to  account  for  the  new  values  of  IPLN. 

Subroutine  UPDATE  was  modified  to  unlock  the  linear  motion  of  a  slip 
joint  based  on  the  parameters  Cl  ar.d  C2  aescriDed  above. 

Subroutine  VISPR  was  modified  to  ac.-onimr.date  the  new  values  of  IPIN. 
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7.0  HYPERELLIPSOID  OPTION 


In  order  to  improve  the  modeling  of  corners  and  other  geometries  the 
option  to  use  hyperellipsoids  as  contact  surfaces  rather  than  standard 
ellipsoids  was  added.  This  option  was  originally  developed  for  General 
Motors  Corporation  and  has  been  incorporated  into  the  ATB  model  with 
their  permission. 


A  hyperellipsoid  is  defined  as  the  surface  generated  by  the  functional: 


F(p)  =  ix/aj®  +  |  y/b| m  +  (z/cj® 


(29) 


where  p 

x,  y,  z 
a,  b,  c 
m 


is  the  vector  from  the  center  of  the  hyperellipsoid  to 

a  point  on  the  boundary, 

are  the  components  of  p, 

are  the  semi-axes'  lengths,  and 

is  the  power  of  the  hyperellipsoid,  an  even  integer. 


As  with  an  ellipsoid  if 


F(p)  -  1  the  point  p  is  on  the  hyperellipsoid  surface, 

F(p)  '  1  the  point  p  is  an  interior  point,  and  if 

F(p)  >  1  the  point  ]p  i6  an  exterior  point. 


if 


If  m  -  2  the  surface  is  an  ellipsoid.  For  larger  values  of  m  the  figure 
squares  off"  at  the  corners.  As  m  approaches  infinity  the  figure 
approaches  a  rectangular  parallelpiped  with  the  same  dimensions  as  the 
hyperellipsoid.  This  makes  the  hyperellipsoid  very  useful  for 
describing  contact  surfaces. 


To  compare  the  hyperellipsoid  shape  to 
rb,  and  z  =  rc  in  the  functional,  then 
point  is  on  the  hyperellipsoid  surface 
powers  of  m. 


a  parallel  piped  let  x  =ra,  y  = 
the  value  for  r  for  which  this 
is  given  in  Table  2  for  various 
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x  =  a,  y  =  b,  and  z  -  c  is  a  vector  to  the  corner  of  the  parallel  piped, 
therefore  r  is  the  ratio  between  the  c cap  orient  b  of  the  vector  on  the 
hy perel 1 i psoi U  surface  in  the  direction  ui  the  parallel  piped  corner  and 
the  components  of  the  vector  co  the  corner- 


Note  that  for  rr  =  128  the  ellipsoid  point  is  within  1/  of  the  corner  of 
the  rectangular  parallelpipod. 


TABLE  2 


Corners  of  a  iiy  perel  1 1  psoid 


M 


r 


6 
lb 
32 
64 
128 
256 
512 
» arge 

7.1  CONTACT  WITH  A  PLANE 


0.5/735 
0.75984 
0 . 37  i  6  9 
0. 93364 
0. 9662.5 
0 . 96296 
0.99145 
0.99572 
0. 997 86 

i-IogiJ)/®  “  1  -  1 . 0 9 86 / m 


Subroutine  PLELP  computes  the  contact  vita  a  p.ane  and  has  been  modified 
to  allow  the  use  of  hyperellipsoids.  If  the  surface  is  an  ellipsoid 
represented  by  the  oid  format  the  original  method  of  calculating  the 
point  of  maximum  penetration  is  used.  If  the  surface  is  a 
hyperel lipsoid  the  surface  point  whose  normal  is  perpendicular  to  the 
plane  is  the  point  of  maximum  penetration  and  can  be  found  by  taking  the 
gradient  of  the  functional  in  equation  29; 
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aT  (30) 

where  a  is  a  positive  scalar  and 

T  is  the  outward  normal  to  the  plane. 

The  coordinates  x,  y  and  z  of  the  point  of  maximum  penetration  are 
readily  obtained  from  equation  30  a6  functions  of  a.  u can  then  be 
computed  by  substituting  these  coordinates  into  equation  29.  This 
computation  is  done  by  the  double  precision  function  HYPEN.  With  the 
value  of  a,  the  point  of  maximum  penetration,  XH  is  computed  and  a  scale 
factor,  FM,  is  determined.  This  scale  factor  time6  the  vector  XH  will 
produce  the  vector  from  the  center  of  the  hyperellipsoid  to  the  plane. 

If  the  surface  is  an  ellipsoid,  thi6  vector  will  locate  the  center  of 
the  intersection  ellipse  in  the  plane.  The  quantity  AMR  =  l-Fff2  is  then 
evaluated  and  if  it  is  less  than  or  equal  to  zero  there  is  no  contact  of 
the  surface  with  the  plane  and  no  further  computations  are  done. 

If  there  is  contact,  the  contact  is  checked  to  determine  if  it  is  within 
the  boundaries  of  the  pLane  and  the  forces  are  applied  as  described  in 
Section  4  of  this  report  and  Vol  I  of  Reference  4. 

It  should  be  noted  that  the  roll-slide  option  can  not  be  used  with 
hy  per  el  1 ipsoi ds . 


7F  =  -  - 


,m-i 


7.2  CONTACT  ANOTHER  H Yl’ERELLIPSOID 


Modifications  were  made  to  subroutine  SEGSEG  to  handle  the  contact  of 
two  hyperellipsoids.  Also  two  new  subroutines  HYEST  and  HYNTR  were 
written  to  replace  subroutine  INTERS  to  calculate  the  penetration 
parameter  and  the  point  of  force  application. 
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To  determine  if  penetration  exists  and  the  point  of  force 
application, three  conditions  must  be  met.  First  the  point  of  force 
application  must  be  the  same  in  both  segment  coordinate  systems: 

q  =  T  (p  -  d).  (31; 


Where 

p  is  the  point  in  the  first  segment's  principal  coordinate 
sy  stan, 

q  is  the  point  in  the  second  segment's  principal  coordinate 
sy  stem, 

d  is  the  vector  between  the  centers  of  the  two  surfaces  in  the 
first  segment's  principal  coordinate  system, 

I  is  the  transformation  matrix  from  the  first  to  the  secona 
segment's  principal  coordinate  6y6tem. 

The  normal  to  each  surtace  passing  through  this  point  must  be  parallel 
and  opposite  in  sign: 


7  F  ( j5 )  =  -  c  .  G  ( q  ) 


(32) 


Where 

F ( p )  is  the  functional  of  the  first  hy perel 1 1 psoi d  at  point  p, 

G(a;  is  the  functional  of  the  second  hy perel 1 i psoid  at  point  a, 

is  the  vector  gradient,  and 

c  is  a  positive  scalar. 

Finally,  the  point  is  chosen  to  be  within  each  hyperellipsoid  a  distance 
proportional  to  the  hyperellipsoid  size: 

F(p)  =  G(q)  (33) 

If  the  value  of  F(p)  (and  G(q))  is  less  than  1  the  figures  intersect,  if 
the  value  is  greater  than  i  no  intersection  occurs  and  if  the  value  is  1 
the  figures  just  touch  at  the  point  p. 
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In  the  original  algorithm,  since  the  figures  are  ellipsoids,  equations 
(29)  and  (30)  may  be  combined  to  form  a  matrix  equation  which  can  be 
solved  for  p  as  a  function  of  the  scalar  c.  A  Newton-Raphson  iterative 
method  is  then  used  to  determine  the  value  of  c  that  allows  all  three 
equations  to  be  satisfied  (Subroutine  INTERS). 

In  the  hy perel 1 ipsoid  case,  equation  30  is  no  longer  a  matrix  equation 
and  so  a  different  approach  was  developed.  To  obtain  a  first 
approximation  the  hyperellipsoids  are  treated  as  rectangular 
paral 1  el  piped  'boxes'  whose  half-widths  are  the  same  as  the  semi-axes  of 
the  hyperellipsoids.  Subroutine  HYEST  determines  whether  or  not  these 
'boxes'  intersect. 

If  the  'boxes'  intersect,  subroutine  HYNTR  is  called  to  refine  the 
estimate . 


7.3  NEW  SUBROUTINES 

A  number  of  new  subroutines  were  added  to  the  ATB  model  in  support  of 
the  hy perel 1 ipsoi d  option.  A  description  of  each  of  these  routines 
follows.  It  should  be  noted  that  in  several  of  these  routines  to  obtain 
approximations  of  desired  quantities  reference  is  made  to  a  'box'.  This 
'box'  is  the  rectangular  paral lelpiped  that  is  centered  at  the  center  of 
a  hyperel 1 ipsoid.  The  edges  are  parallel  to  the  principal  axes  of  the 
hy perel 1 i psoid  and  the  half-widths  are  the  same  as  the  semi-axes  of  the 
hy perel 1 l psoi d .  For  large  values  of  the  power  (the  exponent)  the 
hy perel 1 l psoi d  almost  'fills'  the  box. 

In  order  to  btore  the  variables  required  to  define  a  hyperellipsoid  the 
BD  array  containing  the  ellipsoid  parameter  was  reformatted.  The 
formats  now  used  are  listed  in  Table  3. 
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TABLE  3 


Format  of  BD  Array 


old  format 

new  format 

Subscript 

ellipsoids  only 

hvuerel 1 ipsoi ds 

1 

a 

-1  power  of  x 

2 

h 

a 

3 

c 

b 

A 

0(1) 

c 

5 

0(2) 

0(1) 

b 

0(37 

0(2) 

1 

DtED(1,  i ) 

0(3) 

8 

DtED(2, 1) 

D(l,  1) 

9 

cTED(3, 1 ; 

D(  2,  i  ) 

10 

DTED( 1,2) 

0(3, 1) 

11 

DTEDt  2, 2) 

D(l, 2) 

12 

DtEDU,  1) 

0(2,2) 

13 

DrEDt 1,3) 

D( 3  ,  2 , 

1A 

O1  EDI.  i,3) 

2  U  ,  3  } 

13 

DTEDt3.  3) 

0(2,3) 

16 

DT.-'D(i,  1  ) 

0(3.3) 

17 

DtFDU,  1) 

/  r* 

18 

dtfd( 3 , 1 ; 

Ufcf 

19 

DrFDU,  i) 

1/  C“ 

20 

jUlk  2,  i) 

i  power  of  x 

21 

I^FDt  i.i, 

:n  power  of  y 

22 

D  ■■  F  L ,  i  ,  > ; 

n  power  of  z 

23 

Dv2, 3) 

0  if  equal  pi 

2A 

D-f D(  3 , 3  / 

-used 

where:  i,  m,  n  are  the  power  a  ol  the  hypertl ) lpsoid, 

a,  h,  c  are  the  semi  axes  of  the  ( hy per )el l ipsoi d, 

0  is  the  offset  of  the  ellipsoid  from  the  c.g.  of  tne  segment, 
ii  is  the  direction  cosine  defining  tne  orientation  of  the 

( hy per ) e 1 1 1 psoi c  with  respect  to  the  segment  principal  axes, 
is  the  ellipsoid  matrix,  and 
is  the  inverse  of  the  ellipsoid  matrix. 

7.3.1 _ Subroutine  HYABF  (E,  Z.  A.  i) 


This  routine  computes  the  hy pe rel 1 ipeci d  functional,  F,  and  it's 
derivatives.  It  is  called  by  subroutines  HYEST  and  HYNTR. 


Inputs : 


B  BD  array  for  hy perel 1 i psoi d  containing  m,  a,  b,  c 
Z  array  containing  x,  y,  z 
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Equations : 


F  =  U/a)m  +  (y  /  b)m  +  (z/c)m 

/ym-l/am\ 

(  Vf)/m  =  /  ym-l/ bin 
ym-l/c®/ 

The  diagonal  hyperellipsoid  matrix  is 

(x/aF2  0  0 

0  (y/b)m~2  0 

0  0  (z/c)m_2 

Out  put  s : 

A  3x3  matrix  containing  the  diagonal  elements  of  the 
hyperellipsoid  matrix  in  the  first  column,  the 
components  of  F  in  the  third  column. 

F  the  value  of  the  functional 

7.3.2  S u b rou ti ne _ a YBND( M. Z , IV. U . C. X) 

This  routine  computes  a  point  on  the  polygon,  determined  from  the 
intersection  of  a  plane  with  a  box,  that  is  furthest  from  an  interior 
point  of  the  polygon  in  c  specified  direction.  It  is  called  by 
subroutine  PLEDG. 

Inputs  : 

M  number  of  points  in  array  Z, 

Z  array  determined  by  subroutine  HYBOX, 

IV  pointer  array  determined  by  HYBOX, 


3  7 


U  vector  direction  or  interest , 


C 


r*1' 

L-i, 


use  direction  of  U , 
use  -U  direction, 


Outpuc : 

X  point  on  box  in  direction  of  C*U 
Equat ions : 

The  distance  of  point  j  from  the  origin  in  the  direction 
CU  is  given  b\ 

a  =  cuT2Tvn 


Procedure : 


The  points  are  examined  and  the  one  yielding  the  maximum  d 
is  stored  in  X.  If  two  •  dint  s  give  the  sate  distance  X  is 
their  average  value  i  me  mid- to:  i;t ; . 


7.3. 3 _ subroutine  _  H  YB  OXj,  L  Tj  P ,  Nj  Z1  i  y  j_ 


This  routine  computes  the  intersection  or  a  plane  vith  the  edges  of 
rectangular  box.  It  is  called  hv  subroutine  l-'LEDG. 

Input  s  : 

E  array  containing  a,  b,  c,  the  half-widths  of  the  box. 

T  the  vector  normal,  to  the  plane, 

P  the  coordinates  of  a  point  on  the  plane. 
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Outputs : 


N  the  number  of  paired  solutions, 

IV  pointers  to  on  ordered  set  of  solutions, 

ZC  the  coordinates  of  the  points  of  intersection. 

Equations : 

Let  the  origin  of  coordinates  be  the  center  of  the  box  and  let  the 
vectors  E( i ) ,  i  =  1 ,  2,  3,  be  parallel  to  the  edges  of  the  box  and  of 
length  equal  to  the  respective  half-widths. 

Let  Z  be  the  vector 

Z  -  u  E( 1 )  +  v  E(2)  +  w  E(3),  where  u,  v,  w  are  scalars. 

1  is  a  point  in  the  oox  if  -1  u,  v,  w  1 
It  TTZ  =  T^P  the  point  is  m  the  plane. 


Procedure  : 

The  box  has  b  surfaces,  each  of  these  is  selected  in  turn  and  the 
functional  T^Z  -  T^P  is  evaluated  at  the  four  corners  of  the  surface. 

It  the  functional  changes  sign  between  any  adjacent  corners  the 
plane  intersects  the  edge  of  the  box  between  these  corners. 

The  intersection  point  is  computed  and  stored  in  the  array  ZC. 

for  each  surface,  points  are  obtained  in  pairs,  lines  of 
intersection  ot  zero  length  are  ignored.  The  maximum  number  of  pairs 
will  be  12  and  the  minimum  will  be  b  if  there  is  a  true  intersection  of 
the  plane  with  the  box. 


It  at  least  0  pairs  were  tour.ti  the  location  in  the  array  ZC  of  a 
unique  path  over  the  surface  oi  the  box  is  determined  ana  stored  in  the 
array  IV  in  such  a  fashion  that  the  sequence 

j  =  IV(k)  ,  k  —  i ,  3,  3,...,  N-I  will  determine  the  path  going 
through  the  points  ZC(*, i). 

7.3.4  Subroutine  HYDAD(D, A. DAB; 

This  routine  computes  the  matrix  D^AD  where  D  is  a  direction  cosine 
matrix  and  A  is  a  diagonal  matrix.  It  is  called  by  subroutine  HYNTR. 

Inputs  : 

b  3x3  direction  cosine  matrix, 

A  array  containing  diagonal  elements  of  A  (see  HYABE) 

Output  : 


PAD  the  3x3  product  matrix,  D^Al), 

Procedu re  : 


The  computation  is  a  st  raignt-1.  orwaro  matrix  multiplication. 
h  L-  5. ..  Subroutine _ KYES'KEM,  BH.TAB . 

This  routine  is  caLIed  by  subroutine  SEGSEG  to  make  a  preliminary 
estimation  of  intersection  of  two  hy perel 1 ipsoi ds  if  an  estimate  does 
not  exist.  It  is  called  by  subroutine  SEGSEG. 

Input  : 

BM,  3N  BD  arrays  containing  the  data  defining  the 
hyperel 1 ipsoi ds  m  and  n, 


AO 


The  following  inputs  art  in  COMMON/ TEMPVS/ 


R 

Pi  2 

Output  : 

TAB 

The  toil 

VU. 


the  vector  l  rom  the  center  of  hy perel 1 ipsoid  m  to  the 
center  ot  hy perel 1 l psoi d  n. 

the  direction  cosiue  matrix  which  transforms  for  the 
segment  reference  of  n  to  the  segment  reference  of  m . 


array  used  as  a  memory,  contains  the  same  information 
as  the  V  array  describee  below  if  there  is  an 
i nit r  sect l on . 

owing  outputs  are  in  COMMON /TEMPVS/ 

array  containing  the  tollowing; 

i  value  ot  uvi.P),  ratio  ct  magnitudes  of  gradients  at 
t h e  .  vs terse-i  1 1  o n  p o i  ri t  s , 

»  Value  of  (BE),  tin  expansion  factor, 

. 1  f  H  i  i  I'i  l  0 1 1  iiy  'f.rrtf  ;  i  i  psoi  il  III » 

7  petit  eis  hy  psrellipscid  n. 


hi 


Equat  ioi.s  : 


Let  Z  be  a  vector  on  m  and  U  be  a  vector  on  n  where 

Z  =  U  +  R  ;  i ;;  an  expansion  i  actor, 

,  -  F  ]  /  j  G  |  ,  w.iere  F  and  G  are  the  hy  perel  1  i psoi  d 

functionals  for  Z  and  U  respectively, 

-  Z  -  li  /  ;  R  :  . 


Procedure  : 

The  hy pere 1 l 1 psoi ds  are  treated  as  boxes  and  subroutine  HYLPX  is 
called  to  line  the  largest  value  ot  ,  Z  and  U  that  satisfy  the 
equation  Z  -  U  *  .  P.. 

If  t ne  value  of  is  less  than  1  there  is  no  intersection  and  the 
routine  exits  storing  the  value  ot  .  in  the  TAB  array.  The  routine  also 

exits  for  a  value  of  equal  to  1  since  there  can  be  no  penetration. 

If  the  value  of  is  greater  than  1  there  is  an  intersection  of  the 

boxes.  In  this  case  the  points  Z  end  U  are  scaled  to  lie  on  the 

respective  hy per  el  1 1 psoi ds ,  the  value  ot  ,  u  estimated  and  the  value  of 
for  the  t  cal  ed  points  is  estimated.  The  results  are  stored  in  the  TAB 
a r ray  . 


7.3.  b_  Double  Free i  si on_  Pune t  ion _ HYFC1U  C_,  Z  .  A,  P ) 

This  function  is  used  Dy  subroutines  HYaBF,  HYL1M,  and  HY'VAL  to 
evaluate  the  term  BYFCN  =  C(  Z.  /  A)^  i  u  such  a  fashion  as  to  prevent 
underflows.  The  value  of  A  is  always  greater  than  zero  and  Pis 
non-negative. 


Equations : 


If  P  =  0,  HYFCN  =  C  and  if  Z  =  0,  HYFCN  =  0. 

For  P  >0  and  Z  i  0,  the  vaLue  of  a  parameter  q  is  determined; 


I 

I 

lb 


r 

i 


q  =  P(ln| Z  j  -  InA  ) 

If  q  _  -88.5,  HYFCN  =  0.  Otherwise  HYFCN  =  C  exp(q) 

Procedure  : 


,  The  equations  are  evaluated  as  above.  The  value  of  -88.5  should  be 

J 

I  adjusted  to  represent  the  smallest  value  which  will  not  produce  an 

'  underflow  on  the  computer  being  used.  (exp(-88.5)  =  3.6*10^~39)) 

J  7.3.7  Subroutine _ HYLIM(A.U.B.V.C.W.Z.BD) 

'  This  routine  is  used  to  calculate  the  boundaries  of  the  figure 

J  formed  by  the  intersection  of  a  by  per el  1 ipsoi d  and  a  plane,  i.e.,  the 

.’  point  r„  the  figure  whose  abscissa  is  a  minimum  or  a  maximum.  It  is 

called  by  subroutine  PLEDG. 

Input : 

U  vector  defining  horizontal  axis  (abscissa) 

V  vector  defining  vertical  axis  (ordinate) 

e.  scalar  multiplier  of  Z 

W  vector  such  that  C*VJ  is  in  the  plane 
Z  estimate  ot  desired  point 
BD  array  containing  parameters  of  hy perel 1 l psoid 
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Output 


A  scalar  multiplier  of  l' 

B  scalar  multiplier  ot  V 

Z  desired  point 

Eq  ua 1 1 on  s : 

Z  =  AH  +  BV  CW  ,  origin  is  cencer  ol  hy  per  ellipsoid, 

T  =  OxV  ;  vector  cross  product,  a  vector  perpendicular 
to  plane. 

Oor.st  raint  s  equations, 

2*12  1  ;  ruuctional  equation  ot  h>  perel  1  i  psoi  a, 

VtEZ  -  J  ,  boundary  constraint,  normal  perpendicular  to 
ordinal e, 

T^Z  -•  CT*W  ,  Z  must  lie  ;  r.  plane, 


Procedure : 

A  fust  estimate  was  ottai.net.  before  railing  this  routine  by  e.  call 
to  HYBND.  This  estimate  ts  the  vaiut  of  L  on  entry. 

First  order  perturbation  equations  are-  used  to  refine  the  value  of 
Z.  These  equations  are: 

ZT£D  =  11-  2r£Z)/m 

VT£D  =  -VTfc2/(m-i) 

TTD  =  ctTw  -  T1^ 
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where  D  is  the  perturbation  of  Z  and  m  is  the  power  of  the 
hyperellipsoid. 

The  equations  are  solved  for  D  and  Z  is  updated  as  Z  =  Z  +  D. 


The  process  is  iterated  until  the  perturbations  are  small  when 
compared  to  Z. 

When  convergence  is  obtained  the  values  of  A  and  B  are  computed. 
7.3.8  Subroutine  HYLPR(  J1  J2,  ID.  C,S,E,T) 

Thib  routine  is  a  simplex  method  for  solving  a  linear  programming 
problem.  It  is  called  by  subroutine  HYLPX. 


Input  . 

Ji  index  of  first  column  to  search, 

J2  index  of  last  column  to  6earch, 

ID  pointer  array  to  identify  columns, 

C  cost  vector, 

S  constraint  array, 

E  temporary  storage  for  pivot  col umn, 


Output  : 

T  vector  indicating  final  costs  of  each  column, 
right  hand  column  contains  solutions  obtained, 
pointer  to  identify  columns. 


05 


S 

ID 


Equa  c ions : 


S.  is  a  maLrix  whi.de  rows  are  the  constraint  equations  and 
whose  columns  ate  the  coefficients  of  a  particular 
variable  in  these  constraint  equations.  The  last  column 
ot  S  is  the  constant  term  in  the  equations. 

t(j)  =  ciTs(*, j)  -  cu),  :  -  .11,  J2,  k  =  id(j; 

where  Cl  is  the  cost  vector  of  the  current  solution. 

Procedure : 

The  simplex  algorithm  is  us*--1  The  values  ot  Tlj;  are  computed  if 
any  T(j)  is  positive,  variable  j  is  entered  into  the  solution,  replacing 
the  variable  whose  elimination  '/ill  reduce  the  cost.  Pointers  to  the 
current  solution  variables  are  kept  in  the  ID  array  which  is  updated. 

The  process  is  iterated  until  an  H  j  )  are  non- pos ■ t iv e . 

The  J1  =  J 2  the  variable  identified  with  co'umn  J 1  is  forced  into 
the  solution  ar.u  no  iteialion  is  performed. 

7 .3 .  i Subroutine H  YL  PXt,  b  B  Nj 

This  program  is  called  by  subroutine  tlYEST  to  solve  for  the  estimate 
of  the  points  ot  intersection  of  two  hyperell ipsoids. 

Input  : 

BM,BN  BD  arrays  containing  the  parameters  of  hy perel 1 i psoi ds 


b() 


m  and  n . 


The  following  inputs  are  in  COMMON/ TEMPVS/ 


R  the  vector  from  the  center  of  m  to  the  center  of  n. 


£  the  direction  cosine  matrix  which  transforms  from  the  n's 

reference  system  to  m's  reference  system. 


Output:  (in  COMMON /TEMPVS/ ) 


VU) 


an  array  containing  the  following; 


i  =  1,  2,  3 
i  =  4,  5,  6 
l  =  7 


point  on  the  box  enclosing  hy perel 1 ipsoid  m, 
point  on  the  box  enclosing  hyperellipsoid  n, 
the  expansion  factor. 


Constraint  Equations: 

Z  V  -  .  R  =  U 


ZU). 

a  ( l  ) , 

where  a(i) 

are  the 

semi-axes  ot  m, 

i  =  1,3, 

V(l); 

b(i). 

where  b(i) 

are  the 

semi -axes  ot  n, 

i  =  1,3, 

Procedure : 

The  array  S  representing  the  constraint  equations  is  computed.  The 
value  of  is  assigned  a  cost  of  -1,  the  values  of  Z  and  V  and  their 
associated  slack  vectors  are  assigned  costs  of  0-  Subroutine  HYLPR  is 
called  to  solve  tor  the  values  ot  Z  and  V  which  produce  the  maximum 
value  of  . 

It  atter  the  initial  call  to  HYLPR  the  associated  cost  vector 
indicates  that  there  is  more  than  one  solution  HYLPR  is  recalled  to  find 
all  solutions  and  the  results  are  averaged. 


See  the  descr  1  pt ior.8  of  subroutines  HYEST  and  HYLPR  for  more 
detai L  s . 


Z- 3^  10 _ Subroutine _ HYNTRCBM.BN,  TAB  ) 

This  subroutine  is  called  by  subroutine  SEGSEG  to  determine  the 
points  on  intersecting  hy perel 1 i psoi ds  that  are  used  to  determine  the 
penetration  (if  any)  of  these  figures. 


Inputs : 

BM,  BN  BD  arrays  containing  the  parameters  of  the  figures 

m  and  n, 

TAB  array  containing  the  current  estimates  of  the 

desired  points, 

The  following  inputs  are  in  COMMON /TEMPVS/ 

R  the  vector  from  the  center  of  figure  m  to  the  center 

of  figure  n. 

D12  the  direction  cosine  matrix  which  transforms  a  vector 

from  the  segment  reference  system  associated  with  n 
to  that  of  m, 


Output 


TAB(i)  an  array  containing: 


1=1  the  value  of  i  , 


1=2  the  value  of 


l 


=  .3,  4,  3  the  value  of  Z,  the  point  on  m, 


i  =  6,  7,  8  the  value  of  V,  the  point  on  n. 

This  TAB  array  is  offset  such  that  i  =  1  corresponds  to  a  value  of  i 
=  23  in  the  TAB  array  used  by  subroutine  SEGSEG. 

Equations : 

Z  =  V  +  ,  R;  relation  between  points, 

■  F  =  -  -J.V  G,  where  F  and  G  are  the  hyperellipsoid  functionals 
for  the  m  and  n  hyperellipsoid6  respectively. 

F  =  1,  G  =  1 

Let  F  =  ZTAZ,  G  -  VrB V ,  •  F  =  AZ ,  7 G  =  BV,  and  let  dZ ,  dV,  d  t 

and  di.  be  perturbations  of  2,  V,  ,  and  respectively.  The  linearized 
perturbation  equations  are: 

dZ  -  dV  -  d  .  R  -  -Z  +  V  +  ,  R 

dAZ  +  <  d£V  +  d  ■  fiV  -  -&L  -  .BV 

VTjj  dV  -  i  -  VTjJV 

ZTAdZ  --  1  -  ZT^Z 

Procedure  : 

The  values  stored  m  the  TAB  array  on  entry  are  used  as  first 
guesses  to  the  variables  <  ,  ,  Z  and  V. 

The  perturbation  equations  are  solved  and  the  values  updated.  The 
procedure  is  iterated  until  the  perturbations  of  Z  are  small  compared  to 
the  value  of  Z. 


When  convergence  is  determined  the  TAB  array  is  updated. 
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t 


This  subroutine  is  called  by  subroutine  PLEDG  to  compute  an 
approximate  area  and  centroid  for  the  figure  formed  by  the  intersection 
of  a  hyperellipsoid  and  a  plane. 

All  input o  and  outputs  are  in  COMMON / TEMPVS/ . 

Input  : 

AMI,  AM2,  AFM,  AEP,  AX1,  AX2,  coordinates  of  the  boundaries  of 

BM1,  BM2,  BEM,  BEP,  BX1,  BX2,  the  figure 

AMIN,  AMAX,  RMIN,  BMAX 

Output  : 

AREA  proportional  to  tne  area,  the  true  area  is  this  number 
times  the  magnitude  of  the  cross  product  of  the  vectors 
used  to  detine  the  abscissa  and  the  ordinate  of  the 
coordinate  system  used  in  PL  ETX7. 

AH  the  location  of  the  centroid  m  trie  abscissa  coordinate, 

Bb  tiu  location  ot  the  centroid  in  the  ordinate  coordinate, 

Ec  uat  l  uiis  : 

Consider  the  area  below  the  straight  iiae  segment  for  the  point 
(xl,yl)  to  the  point  (x2,y2).  Then 

dx  =  x2  -  xl 

ar  =  dx(y2  +  yl),  twice  the  increment  of  area, 

ax  =  ar(.x2  +  xil  +  dx(x2y2  +  xlyl),  six  times  the  increment 
of  the  abscissa  of  the  centroid, 
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ay  =  ar(y2  +  yl)  -  dxy2yl,  six  limes  the  increment  of  the 
ordinate  of  the  centroid. 


AREA  =  sum  of  the  ar  divided  by  two, 


sum  of 

the 

ax 

div l ded 

by 

six  times 

the  area 

sum  ol 

the 

ay 

divided 

by 

six  times 

the  area 

Procedure : 

The  sections  ot  the  general  shape  are  shown  in  Figure  y. 

Tests  are  made  lor  the  existence  of  the  sections  and  the  formulae  are 
used  to  compute  the  area  and  the  centroid. 


(A” : 
<  RM  i 


(ARM ,  BMAX  ) 
(A K..'  .  UMAX)  i 


7 


/ 

/ 


/ 


:na:-:  (  i, A  !  ,  BM  I  A ) ) 
mint  HA 1 , HMAX ) ) 


/ 

\ 


\ 


\ 

( HM2 , BM , M  • 
l  HKM.HMIX  ) 


«  (  ARP  ,  BMA.X  ; 


\ 


(AX! , BMAX ) 


\ 

\ 


\ 

/ 


( AX2 , max ( BA2 , BM1 N ) > 
(BX2 , m i n ( BA2 , BMAX ) ) 


/ 


/ 


/  (BX.1  ,  HM  IN) 
-•  (HFP.BMIN) 


{ .lb.’. i  i  as. i ,  ■  ■  r d  i  11.1 1  r  ) 


Fi  i'll  re  ')  Hyperel  I  ijisoid  Romrion  Area  Boundaries 
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_7_-  3 2_  _  Subrouc  1  lit  H  VSUL  (, A,  N ,  N  i)) 


This  routine  is  a  modiiieti  Gauss  Eiimindtion  process  tor  solving  a 
set  or  simultaneous  equations.  It  is  called  by  subroutine  HYNTR. 


Input  : 


A  array  containing  the  simultaneous  equations, 


N  thf  number  of  equations, 


NT  the  tirst  subscript  01  the  array  A.. 


Output  : 


A.  the  reduced  (solved  equal  tons.  The  solutions  art  in  column 
Nei  ot  the  array. 


Proce  dure  : 


Gauss  -  Elimination  ;s  used  with  the  pivot  always  cn  the  diagonal. 
No  pivoting  1  .■>  done  for  a  z<  i  c  d.ag.nal  .  This  modilication  was 
necessary  beeiiUu ,  the  m  a :  \  1  x  may  nave  at.  a  .  i  r  c  r  o  c  o  1  tun  n  (with  a 
corresponding  al;  zerc  row  . 

7.3.13 _ Subroutine.  K  YVAl.(  Ajl'j  R,3  Jr,  I. ; 

This  routine  ccmpu'es  the  point  on  a  hvporel  ilpsoi  c  that  lies  on  a 
particular  Line.  1 1  is  called  by  subroutine  PL  EDO . 

Input  : 

U  vector  defining  line  ot  interest, 

R  vector  locating  end  point  of  line, 
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BD  array  containing  hyperel  1  ipsoid  parameters, 

L  =  1  indicates  point  along  -U  desired, 

=  2  indicates  point  along  +0  desired, 

Output : 

A  parameter  detimng  point  on  U. 

£q  ua  t ions : 

Z  -  AU  +  R,  point  on  line  U, 

zTbz  1,  constraint  that  point  lies  on  hyperellipsoid, 

1 hy pore  1 1 1 psoid  iunctional) 

t  t  ui  bat  l  on  equat  j  ons  : 

Let  I  l  A ,»  --  ZTfiZ  -  1,  and  let  e  be  a  perturbation  ot  A.  Expand  up 

to  the  second  derivative  with  respect  to  A  to  get 

I(A+e)  =■  tl.A)  +  el  '  +  e~t "/  2  =  0, 

VI'  is  the  tirst  derivative  ot  t,  and  t"  is  the  seconc  at  A), 
e  -  -i  (  A )  /  1 1  '  /  2  e  s  i  gni  £  '  )  l  U  '  /  2  ) 2  -  f  (  A )  f  "/ 2  J 1  /  2  ] 


Procedure : 


Subroutine  H YVB A  is  called  to  determine  the  point  on  the  box 
surrounding  the  hyperel 1 l psoid  in  the  direction  specified  by  L. 

Ihe  quadratic  perturbation  equations  are  solveu  and  the  process 
iterated  until  the  tunctional  equation  satisfies  a  prescribed  test. 


I  1 


7_. _3_.J.4_  _  Subroutine _ HYVBX(Q  ,  S .  B  ,  H,  RM/ 

This  routine  is  called  by  subroutine  HYVAL  to  determine  the  points 
ot  intersection  ot  a  vector  with  a  box. 

Input  : 

Q  vector  which  intersects  box, 

S  iixrd  vector  from  c e n t < r  of  box  to  Q 
P  array  containing,  iinu-nsums  o!  box, 

Output  ■ 

M  index  identifying  solutions,  M  -  Z  it  there  is  an 
intersection,  M  -  u  tot  none. 

K>!  array  emit  a:  a  mg  ;  tie  i«\,  poi  tu  »  oi  intersection  if  they  exist, 
fc.quie  ions: 

'L  •••  Af  +  S,  i.eueral  >i  l  nt  on  •> , 


P ri cedu  re : 


It  Z  is  a  point  on  a  face  ot  the  box  some  component  ot.  Z  must  equal 
tnc  dimension  ot  Liie  box.  •:  •  t 


hli;  be  ■»  or  -  trie  halt  widths  ot  the  box, 


Z ( l ) ,  Q ( l )  and  Sf i  )  be  the  corresponding,  components  ot  the 

vector  s  X ,  Q  , 


and  b  respectively. 


11  ,  ibU) 


S 1 1 ;  v  t  )  f 


s',  j  m  1 1 J, 


B 1  j  )Q  ( l  ),  t  hen 


the  point  determined  by  r  -  i.B(i  ;  -  S(i)  )/Q(i)  will  be  a  desired 
point.  The  inequality  is  tested  tor  all  combinations  of  the  indites  and 
the  unique  solution.!;  tot  r  are  saved  in  the  array  RM.  Before  exiting 
the  RM  1  s  are  ordered  such  that  RM(.  1  )  ts  the  smallest  algebraic  solution 
and  RM‘,  1)  is  the  largest. 


7  .  3  .  1  ’j  tunc  tun 


HYRENUDj  r,  V j 


This  routine  js  ,  ailed  by  subroutine  PL.E1.R  to  compute  the  value  ol 
Al.F  which  is  used  m  the  computation  oi  the  intersection  ut  a 
by  pert  1  i  i  psoi  c.  with  a  plane.  The  powers  of  the  hyperel  lipsoid 

1  tiiu'i  !  ri'u  i  r.’t;  i;;  1  i  . 


.2  :  •.■■lit  ,1  :  n  i  L)  r  t  r<  l  i  i  psoi  c.  ■  r.t  uncat  1  on, 


.ir.  d  l  .in, 


v. OV: . i, i, t  d non. 


i  ,  tin-  h  y  p  t  ■  r  (  1  I  j  ;  m  •  1 

nib-  l.i  t  f  I  V  d  t.  i  i  !1  #  (  ■  V 


^  l->  I  ’i/i  t  K  , 


oil  (J 


At  the  poi  iit  ot  maximum  penetration  the  gradient  must  be  parallel  to 
the  plane  vector  1  t  p  ,  t  ,  1 3  )  ,  i.c. 

x  1  a .  K  ^  ~  Ae.  r  L  |  a  /  k  , 

,  y /  b .  '  ~ =•  ALP,  .  b.<  1  . 

,  z/c  r,)"i  ALP.  t_s.  C/K. . 

trt  ■  v>  J  -  1  t  j  a,’  k  ]  l'l  ,  e  ,  -  K/  (  k-  1  /  , 

V  _  '  !  ,t:,  L'/  1  !  lV,  e_.  -  i  ■  ;  ■■  ’  , 

wrier.  •  i.  ■■  ire  computed  by 


:  ,n'.  ,  ■  :  :  stct!  o.  .  c-  ci.ii.tM  i  ba  ‘..lie  lunetionai 

•q  tuit  i«‘i  ■  t  r  t  ,  1  r  s  1  -  ca  r  1  :  ..  TK 1  >.  equation  may  bt: 

n  1  • 

i'  1  >J  b  {'  ■  :  ’  ■■  f.J .  f  '  *  V  ■  It*  .  . 

p  t  oCt-tlu  r  e  . 


t  n  e  n:  a  c  i  n.  uni 


1  lit-  i  xponi  nt  , 


The  value  of  F(ALP)  is  computed.  If  lF(ALPJj  is  less  than  10-^  the 
functional  is  assumed  to  be  satisfied  and  the  routine  exits.  If  the 
t'(ALP)  is  greater  or  equal  to  10“®  a  stepping  procedure  is  used  to 
modify  ALP  until  convergence  is  obtained.  Note:  if  the  exponents  are 
all  equal,  the  first  estimate  of  ALP  should  satisfy  the  convergence 
test. 


8.0  OTHER  NEW  OPTIONS 


A  number  of  minor  modifications  and  corrections  have  been  made  to  the 
ATB1II  model.  These  changes  have  a  label  other  than  the  subroutine  name 
in  column  73  of  the  program  listing  in  Volume  3.  Most  of  these  changes 
have  been  to  correct  errors  in  the  code,  to  add  stops  to  avoid  input 
errors  or  to  improve  the  output  format.  Those  changes  that  allow  the 
user  to  choose  a  new  option  are  described  briefly  below,  A  more 
complete  description  of  the  use  of  these  options  1 s  in  the  input 
description  in  Volume  2. 

Three  new  time  histories  have  been  added.  The  H.8  cards  allow  the  wind 
lurces  on  any  segment  to  be  output,  the  total  forces  and  torques  at  a 
joint  can  be  output  using  the  H.9  cards,  and  body  properties  of  a  single 
segment  or  a  set  of  segments  can  now  be  output  using  the  H.10  cards. 
These  body  properties  include  the  center  ot  mass  location,  total  linear 
and  angular  momentum  and  kinetic  energy.  The  H.l  cards  now  allow  the 
acceleration  output  to  include  the  effects  of  gravity.  This  allows  the 
accelerations  to  be  compared  exactly  with  accelerometer  data.  Ai  ."o  Liu 
user  now  has  ttie  option  to  choose  the  reference  system  in  which  many  ot 
tne  time  hi  slot les  are  output,  by  specifying  KREF  m  the  H  cards. 
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